{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OSL2CMEzmQvB"
   },
   "source": [
    "# Homework 2 - Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ox7joE3aZkh-"
   },
   "source": [
    "Binary classification is one of the most fundamental problem in machine learning. In this tutorial, you are going to build linear binary classifiers to predict whether the income of an indivisual exceeds 50,000 or not. We presented a discriminative and a generative approaches, the logistic regression(LR) and the linear discriminant anaysis(LDA). You are encouraged to compare the differences between the two, or explore more methodologies. Although you can finish this tutorial by simpliy copying and pasting the codes, we strongly recommend you to understand the mathematical formulation first to get more insight into the two algorithms. Please find [here](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Logistic%20Regression%20(v3).pdf) and [here](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Classification%20(v3).pdf) for more detailed information about the two algorithms.\n",
    "\n",
    "二元分類是機器學習中最基礎的問題之一，在這份教學中，你將學會如何實作一個線性二元分類器，來根據人們的個人資料，判斷其年收入是否高於 50,000 美元。我們將以兩種方法: logistic regression 與 generative model，來達成以上目的，你可以嘗試了解、分析兩者的設計理念及差別。針對這兩個演算法的理論基礎，可以參考李宏毅老師的教學投影片 [logistic regression](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Logistic%20Regression%20(v3).pdf) 與 [generative model](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Classification%20(v3).pdf)。\n",
    "\n",
    "若有任何問題，歡迎來信至助教信箱 ntu-ml-2020spring-ta@googlegroups.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nkNW5cQmohoo"
   },
   "source": [
    "# Dataset\n",
    "\n",
    "This dataset is obtained by removing unnecessary attributes and balancing the ratio between positively and negatively labeled data in the [**Census-Income (KDD) Data Set**](https://archive.ics.uci.edu/ml/datasets/Census-Income+(KDD)), which can be found in [**UCI Machine Learning Repository**](https://archive.ics.uci.edu/ml/index.php). Only preprocessed and one-hot encoded data (i.e. *X_train*,  *Y_train* and *X_test*) will be used in this tutorial. Raw data (i.e. *train.csv* and *test.csv*) are provided to you in case you are interested in it.\n",
    "\n",
    "這個資料集是由 [**UCI Machine Learning Repository**](https://archive.ics.uci.edu/ml/index.php) 的 [**Census-Income (KDD) Data Set**](https://archive.ics.uci.edu/ml/datasets/Census-Income+(KDD)) 經過一些處理而得來。為了方便訓練，我們移除了一些不必要的資訊，並且稍微平衡了正負兩種標記的比例。事實上在訓練過程中，只有 X_train、Y_train 和 X_test 這三個經過處理的檔案會被使用到，train.csv 和 test.csv 這兩個原始資料檔則可以提供你一些額外的資訊。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 221
    },
    "colab_type": "code",
    "id": "Ww4-VJoJqE-_",
    "outputId": "b969b581-fd83-4bfd-fccd-3b2991b20811"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/bin/sh: 1: gdown: not found\n",
      "data/\n",
      "data/sample_submission.csv\n",
      "data/test_no_label.csv\n",
      "data/train.csv\n",
      "data/X_test\n",
      "data/X_train\n",
      "data/Y_train\n",
      "acc.png\t\t\t       hw2_classification.ipynb  output_generative.csv\n",
      "data\t\t\t       hw2_classification.py\t output_logistic.csv\n",
      "data.tar.gz\t\t       hw2_slides.pptx\n",
      "hw2_classification_demo.ipynb  loss.png\n"
     ]
    }
   ],
   "source": [
    "!gdown --id '1KSFIRh0-_Vr7SdiSCZP1ItV7bXPxMD92' --output data.tar.gz\n",
    "!tar -zxvf data.tar.gz\n",
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WRXI0kf0W4Bd"
   },
   "source": [
    "# Logistic Regression\n",
    "\n",
    "In this section we will introduce logistic regression first. We only present how to implement it here, while mathematical formulation and analysis will be omitted. You can find more theoretical detail in [Prof. Lee's lecture](https://www.youtube.com/watch?v=hSXFuypLukA).\n",
    "\n",
    "首先我們會實作 logistic regression，針對理論細節說明請參考[李宏毅老師的教學影片](https://www.youtube.com/watch?v=hSXFuypLukA)\n",
    "\n",
    "### Preparing Data\n",
    "\n",
    "Load and normalize data, and then split training data into training set and development set.\n",
    "\n",
    "下載資料，並且對每個屬性做正規化，處理過後再將其切分為訓練集與發展集。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "colab_type": "code",
    "id": "7NzAmkzU2MAS",
    "outputId": "07d466eb-8b83-47d8-894a-0d990407afb3"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of training set: 48830\n",
      "Size of development set: 5426\n",
      "Size of testing set: 27622\n",
      "Dimension of data: 510\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.random.seed(0)\n",
    "X_train_fpath = './data/X_train'\n",
    "Y_train_fpath = './data/Y_train'\n",
    "X_test_fpath = './data/X_test'\n",
    "output_fpath = './output_{}.csv'\n",
    "\n",
    "# Parse csv files to numpy array\n",
    "with open(X_train_fpath) as f:\n",
    "    next(f)\n",
    "    X_train = np.array([line.strip('\\n').split(',')[1:] for line in f], dtype = float)\n",
    "with open(Y_train_fpath) as f:\n",
    "    next(f)\n",
    "    Y_train = np.array([line.strip('\\n').split(',')[1] for line in f], dtype = float)\n",
    "with open(X_test_fpath) as f:\n",
    "    next(f)\n",
    "    X_test = np.array([line.strip('\\n').split(',')[1:] for line in f], dtype = float)\n",
    "\n",
    "def _normalize(X, train = True, specified_column = None, X_mean = None, X_std = None):\n",
    "    # This function normalizes specific columns of X.\n",
    "    # The mean and standard variance of training data will be reused when processing testing data.\n",
    "    #\n",
    "    # Arguments:\n",
    "    #     X: data to be processed\n",
    "    #     train: 'True' when processing training data, 'False' for testing data\n",
    "    #     specific_column: indexes of the columns that will be normalized. If 'None', all columns\n",
    "    #         will be normalized.\n",
    "    #     X_mean: mean value of training data, used when train = 'False'\n",
    "    #     X_std: standard deviation of training data, used when train = 'False'\n",
    "    # Outputs:\n",
    "    #     X: normalized data\n",
    "    #     X_mean: computed mean value of training data\n",
    "    #     X_std: computed standard deviation of training data\n",
    "\n",
    "    if specified_column == None:\n",
    "        specified_column = np.arange(X.shape[1])\n",
    "    if train:\n",
    "        X_mean = np.mean(X[:, specified_column] ,0).reshape(1, -1)\n",
    "        X_std  = np.std(X[:, specified_column], 0).reshape(1, -1)\n",
    "\n",
    "    X[:,specified_column] = (X[:, specified_column] - X_mean) / (X_std + 1e-8)\n",
    "     \n",
    "    return X, X_mean, X_std\n",
    "\n",
    "def _train_dev_split(X, Y, dev_ratio = 0.25):\n",
    "    # This function spilts data into training set and development set.\n",
    "    train_size = int(len(X) * (1 - dev_ratio))\n",
    "    return X[:train_size], Y[:train_size], X[train_size:], Y[train_size:]\n",
    "\n",
    "# Normalize training and testing data\n",
    "X_train, X_mean, X_std = _normalize(X_train, train = True)\n",
    "X_test, _, _= _normalize(X_test, train = False, specified_column = None, X_mean = X_mean, X_std = X_std)\n",
    "    \n",
    "# Split data into training set and development set\n",
    "dev_ratio = 0.1\n",
    "X_train, Y_train, X_dev, Y_dev = _train_dev_split(X_train, Y_train, dev_ratio = dev_ratio)\n",
    "\n",
    "train_size = X_train.shape[0]\n",
    "dev_size = X_dev.shape[0]\n",
    "test_size = X_test.shape[0]\n",
    "data_dim = X_train.shape[1]\n",
    "print('Size of training set: {}'.format(train_size))\n",
    "print('Size of development set: {}'.format(dev_size))\n",
    "print('Size of testing set: {}'.format(test_size))\n",
    "print('Dimension of data: {}'.format(data_dim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "imgCeBDoApdb"
   },
   "source": [
    "\n",
    "### Some Useful Functions\n",
    "\n",
    "Some functions that will be repeatedly used when iteratively updating the parameters.\n",
    "\n",
    "這幾個函數可能會在訓練迴圈中被重複使用到。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hSDAw5LTAs2o"
   },
   "outputs": [],
   "source": [
    "def _shuffle(X, Y):\n",
    "    # This function shuffles two equal-length list/array, X and Y, together.\n",
    "    randomize = np.arange(len(X))\n",
    "    np.random.shuffle(randomize)\n",
    "    return (X[randomize], Y[randomize])\n",
    "\n",
    "def _sigmoid(z):\n",
    "    # Sigmoid function can be used to calculate probability.\n",
    "    # To avoid overflow, minimum/maximum output value is set.\n",
    "    return np.clip(1 / (1.0 + np.exp(-z)), 1e-8, 1 - (1e-8))\n",
    "\n",
    "def _f(X, w, b):\n",
    "    # This is the logistic regression function, parameterized by w and b\n",
    "    #\n",
    "    # Arguements:\n",
    "    #     X: input data, shape = [batch_size, data_dimension]\n",
    "    #     w: weight vector, shape = [data_dimension, ]\n",
    "    #     b: bias, scalar\n",
    "    # Output:\n",
    "    #     predicted probability of each row of X being positively labeled, shape = [batch_size, ]\n",
    "    return _sigmoid(np.matmul(X, w) + b)\n",
    "\n",
    "def _predict(X, w, b):\n",
    "    # This function returns a truth value prediction for each row of X \n",
    "    # by rounding the result of logistic regression function.\n",
    "    return np.round(_f(X, w, b)).astype(np.int)\n",
    "    \n",
    "def _accuracy(Y_pred, Y_label):\n",
    "    # This function calculates prediction accuracy\n",
    "    acc = 1 - np.mean(np.abs(Y_pred - Y_label))\n",
    "    return acc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OxJdfhEEOYwg"
   },
   "source": [
    "### Functions about gradient and loss\n",
    "\n",
    "Please refers to [Prof. Lee's lecture slides](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Logistic%20Regression%20(v3).pdf)(p.12) for the formula of gradient and loss computation.\n",
    "\n",
    "請參考[李宏毅老師上課投影片](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Logistic%20Regression%20(v3).pdf)第 12 頁的梯度及損失函數計算公式。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DqYkUgLjOWi1"
   },
   "outputs": [],
   "source": [
    "def _cross_entropy_loss(y_pred, Y_label):\n",
    "    # This function computes the cross entropy.\n",
    "    #\n",
    "    # Arguements:\n",
    "    #     y_pred: probabilistic predictions, float vector\n",
    "    #     Y_label: ground truth labels, bool vector\n",
    "    # Output:\n",
    "    #     cross entropy, scalar\n",
    "    cross_entropy = -np.dot(Y_label, np.log(y_pred)) - np.dot((1 - Y_label), np.log(1 - y_pred))\n",
    "    return cross_entropy\n",
    "\n",
    "def _gradient(X, Y_label, w, b):\n",
    "    # This function computes the gradient of cross entropy loss with respect to weight w and bias b.\n",
    "    y_pred = _f(X, w, b)\n",
    "    pred_error = Y_label - y_pred\n",
    "    w_grad = -np.sum(pred_error * X.T, 1)\n",
    "    b_grad = -np.sum(pred_error)\n",
    "    return w_grad, b_grad\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XXEFuqydaA34"
   },
   "source": [
    "### Training\n",
    "\n",
    "Everything is prepared, let's start training! \n",
    "\n",
    "Mini-batch gradient descent is used here, in which training data are split into several mini-batches and each batch is fed into the model sequentially for losses and gradients computation. Weights and bias are updated on a mini-batch basis.\n",
    "\n",
    "Once we have gone through the whole training set,  the data have to be re-shuffled and mini-batch gradient desent has to be run on it again. We repeat such process until max number of iterations is reached.\n",
    "\n",
    "我們使用小批次梯度下降法來訓練。訓練資料被分為許多小批次，針對每一個小批次，我們分別計算其梯度以及損失，並根據該批次來更新模型的參數。當一次迴圈完成，也就是整個訓練集的所有小批次都被使用過一次以後，我們將所有訓練資料打散並且重新分成新的小批次，進行下一個迴圈，直到事先設定的迴圈數量達成為止。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 255
    },
    "colab_type": "code",
    "id": "s6yNUeG9aBR1",
    "outputId": "5a82e26f-f387-4b89-f15e-1bd45d4f8521"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_train_pred1:\n",
      " [5.29192971e-03 5.89540747e-01 1.00000000e-08 ... 1.71167407e-03\n",
      " 3.41123771e-03 3.45270473e-02]\n",
      "y_pred:\n",
      " [5.29192971e-03 5.89540747e-01 1.00000000e-08 ... 1.71167407e-03\n",
      " 3.41123771e-03 3.45270473e-02]\n",
      "lb:\n",
      " [ -5.24157232  -0.52841144 -18.42068074 ...  -6.3702834   -5.68068009\n",
      "  -3.36601228]\n",
      "y_train_pred2:\n",
      " [5.29192971e-03 5.89540747e-01 1.00000000e-08 ... 1.71167407e-03\n",
      " 3.41123771e-03 3.45270473e-02]\n",
      "y_pred:\n",
      " [1.71636689e-01 9.30497876e-03 1.90581337e-01 ... 1.00000000e-08\n",
      " 3.06454312e-02 3.10066970e-03]\n",
      "lb:\n",
      " [ -1.76237531  -4.67720567  -1.65767621 ... -18.42068074  -3.48527169\n",
      "  -5.77613716]\n",
      "y_train_pred1:\n",
      " [3.75322694e-02 8.04661919e-04 1.00000000e-08 ... 1.36930599e-02\n",
      " 2.05959472e-01 3.18925531e-03]\n",
      "y_pred:\n",
      " [3.75322694e-02 8.04661919e-04 1.00000000e-08 ... 1.36930599e-02\n",
      " 2.05959472e-01 3.18925531e-03]\n",
      "lb:\n",
      " [ -3.2825542   -7.12508835 -18.42068074 ...  -4.29086615  -1.58007587\n",
      "  -5.74796784]\n",
      "y_train_pred2:\n",
      " [3.75322694e-02 8.04661919e-04 1.00000000e-08 ... 1.36930599e-02\n",
      " 2.05959472e-01 3.18925531e-03]\n",
      "y_pred:\n",
      " [2.40592857e-01 1.30463138e-02 1.44876711e-01 ... 1.00000000e-08\n",
      " 1.48082883e-01 8.35924830e-03]\n",
      "lb:\n",
      " [ -1.42464916  -4.33924965  -1.93187216 ... -18.42068074  -1.90998314\n",
      "  -4.78438677]\n",
      "y_train_pred1:\n",
      " [2.00186815e-01 6.59202288e-07 8.11639686e-01 ... 9.74761754e-01\n",
      " 2.73720453e-04 3.57397924e-01]\n",
      "y_pred:\n",
      " [2.00186815e-01 6.59202288e-07 8.11639686e-01 ... 9.74761754e-01\n",
      " 2.73720453e-04 3.57397924e-01]\n",
      "lb:\n",
      " [ -1.60850428 -14.23223539  -0.20869877 ...  -0.02556219  -8.20340322\n",
      "  -1.02890548]\n",
      "y_train_pred2:\n",
      " [2.00186815e-01 6.59202288e-07 8.11639686e-01 ... 9.74761754e-01\n",
      " 2.73720453e-04 3.57397924e-01]\n",
      "y_pred:\n",
      " [1.81709543e-01 1.87821734e-02 1.81930002e-01 ... 1.00000000e-08\n",
      " 4.03464150e-02 7.76014372e-03]\n",
      "lb:\n",
      " [ -1.70534578  -3.97484708  -1.70413327 ... -18.42068074  -3.21025273\n",
      "  -4.85875442]\n",
      "y_train_pred1:\n",
      " [0.01110497 0.81434288 0.51730305 ... 0.01459968 0.07030657 0.02770799]\n",
      "y_pred:\n",
      " [0.01110497 0.81434288 0.51730305 ... 0.01459968 0.07030657 0.02770799]\n",
      "lb:\n",
      " [-4.50036227 -0.20537377 -0.65912641 ... -4.22675582 -2.65489005\n",
      " -3.58603458]\n",
      "y_train_pred2:\n",
      " [0.01110497 0.81434288 0.51730305 ... 0.01459968 0.07030657 0.02770799]\n",
      "y_pred:\n",
      " [2.22121617e-01 1.16160715e-02 1.46374579e-01 ... 1.00000000e-08\n",
      " 1.01117731e-01 9.95612529e-03]\n",
      "lb:\n",
      " [ -1.50453022  -4.45536566  -1.92158633 ... -18.42068074  -2.29146978\n",
      "  -4.60956731]\n",
      "y_train_pred1:\n",
      " [1.49353566e-02 3.79087719e-01 9.48973968e-02 ... 1.25968048e-01\n",
      " 3.27600067e-04 1.24888162e-02]\n",
      "y_pred:\n",
      " [1.49353566e-02 3.79087719e-01 9.48973968e-02 ... 1.25968048e-01\n",
      " 3.27600067e-04 1.24888162e-02]\n",
      "lb:\n",
      " [-4.20402395 -0.96998765 -2.35495901 ... -2.071727   -8.023717\n",
      " -4.38292174]\n",
      "y_train_pred2:\n",
      " [1.49353566e-02 3.79087719e-01 9.48973968e-02 ... 1.25968048e-01\n",
      " 3.27600067e-04 1.24888162e-02]\n",
      "y_pred:\n",
      " [1.97365713e-01 1.89451840e-02 1.56763872e-01 ... 1.00000000e-08\n",
      " 7.71980956e-02 1.17798067e-02]\n",
      "lb:\n",
      " [ -1.62269686  -3.96620552  -1.85301461 ... -18.42068074  -2.56138049\n",
      "  -4.44136851]\n",
      "y_train_pred1:\n",
      " [1.48961290e-02 3.32036836e-01 4.77236612e-04 ... 4.10325922e-01\n",
      " 6.52418129e-02 1.00000000e-08]\n",
      "y_pred:\n",
      " [1.48961290e-02 3.32036836e-01 4.77236612e-04 ... 4.10325922e-01\n",
      " 6.52418129e-02 1.00000000e-08]\n",
      "lb:\n",
      " [ -4.2066539   -1.10250936  -7.64749815 ...  -0.8908035   -2.72965471\n",
      " -18.42068074]\n",
      "y_train_pred2:\n",
      " [1.48961290e-02 3.32036836e-01 4.77236612e-04 ... 4.10325922e-01\n",
      " 6.52418129e-02 1.00000000e-08]\n",
      "y_pred:\n",
      " [1.69433959e-01 3.78177126e-02 2.14037119e-01 ... 1.00000000e-08\n",
      " 8.37954775e-02 1.18548091e-02]\n",
      "lb:\n",
      " [ -1.77529205  -3.2749777   -1.54160583 ... -18.42068074  -2.47937624\n",
      "  -4.43502166]\n",
      "y_train_pred1:\n",
      " [0.81951074 0.0128588  0.08968052 ... 0.50236077 0.14871771 0.88226945]\n",
      "y_pred:\n",
      " [0.81951074 0.0128588  0.08968052 ... 0.50236077 0.14871771 0.88226945]\n",
      "lb:\n",
      " [-0.19904777 -4.35372687 -2.41150166 ... -0.68843675 -1.90570535\n",
      " -0.12525777]\n",
      "y_train_pred2:\n",
      " [0.81951074 0.0128588  0.08968052 ... 0.50236077 0.14871771 0.88226945]\n",
      "y_pred:\n",
      " [2.20099639e-01 1.81082630e-02 1.80099204e-01 ... 1.00000000e-08\n",
      " 9.94070376e-02 8.27464822e-03]\n",
      "lb:\n",
      " [ -1.51367493  -4.01138692  -1.71424745 ... -18.42068074  -2.30853237\n",
      "  -4.79455887]\n",
      "y_train_pred1:\n",
      " [5.88070670e-01 9.25378790e-01 2.27422874e-01 ... 5.08447882e-02\n",
      " 2.65592903e-02 7.45333983e-05]\n",
      "y_pred:\n",
      " [5.88070670e-01 9.25378790e-01 2.27422874e-01 ... 5.08447882e-02\n",
      " 2.65592903e-02 7.45333983e-05]\n",
      "lb:\n",
      " [-0.53090815 -0.07755212 -1.48094411 ... -2.97897766 -3.62837567\n",
      " -9.50426323]\n",
      "y_train_pred2:\n",
      " [5.88070670e-01 9.25378790e-01 2.27422874e-01 ... 5.08447882e-02\n",
      " 2.65592903e-02 7.45333983e-05]\n",
      "y_pred:\n",
      " [2.37633356e-01 3.50708402e-02 2.11799461e-01 ... 1.00000000e-08\n",
      " 8.57585486e-02 1.06355788e-02]\n",
      "lb:\n",
      " [ -1.43702631  -3.35038526  -1.55211539 ... -18.42068074  -2.45621951\n",
      "  -4.54355041]\n",
      "y_train_pred1:\n",
      " [0.1596823  0.05834685 0.0486895  ... 0.91108943 0.10730049 0.00134032]\n",
      "y_pred:\n",
      " [0.1596823  0.05834685 0.0486895  ... 0.91108943 0.10730049 0.00134032]\n",
      "lb:\n",
      " [-1.83456905 -2.84134987 -3.02229184 ... -0.09311421 -2.23212207\n",
      " -6.61484936]\n",
      "y_train_pred2:\n",
      " [0.1596823  0.05834685 0.0486895  ... 0.91108943 0.10730049 0.00134032]\n",
      "y_pred:\n",
      " [3.36316871e-01 2.59545876e-02 1.05184548e-01 ... 1.00000000e-08\n",
      " 5.94820233e-02 1.16471081e-02]\n",
      "lb:\n",
      " [ -1.08970149  -3.6514069   -2.25203887 ... -18.42068074  -2.82208114\n",
      "  -4.45269736]\n",
      "y_train_pred1:\n",
      " [4.03589345e-05 4.12341010e-01 9.63016077e-02 ... 1.19569202e-01\n",
      " 3.29485875e-02 5.33740842e-01]\n",
      "y_pred:\n",
      " [4.03589345e-05 4.12341010e-01 9.63016077e-02 ... 1.19569202e-01\n",
      " 3.29485875e-02 5.33740842e-01]\n",
      "lb:\n",
      " [-10.11769776  -0.88590458  -2.34027027 ...  -2.12385998  -3.41280689\n",
      "  -0.62784487]\n",
      "y_train_pred2:\n",
      " [4.03589345e-05 4.12341010e-01 9.63016077e-02 ... 1.19569202e-01\n",
      " 3.29485875e-02 5.33740842e-01]\n",
      "y_pred:\n",
      " [2.40923241e-01 2.06517391e-02 1.69607954e-01 ... 1.00000000e-08\n",
      " 4.63274909e-02 1.07384422e-02]\n",
      "lb:\n",
      " [ -1.4232769   -3.87995575  -1.77426566 ... -18.42068074  -3.07201974\n",
      "  -4.53392524]\n",
      "Training loss: 0.271355435246406\n",
      "Development loss: 0.2896359675026287\n",
      "Training accuracy: 0.8836166291214418\n",
      "Development accuracy: 0.8733873940287504\n"
     ]
    }
   ],
   "source": [
    "# Zero initialization for weights ans bias\n",
    "w = np.zeros((data_dim,)) \n",
    "b = np.zeros((1,))\n",
    "\n",
    "# Some parameters for training    \n",
    "max_iter = 10\n",
    "batch_size = 8\n",
    "learning_rate = 0.2\n",
    "\n",
    "# Keep the loss and accuracy at every iteration for plotting\n",
    "train_loss = []\n",
    "dev_loss = []\n",
    "train_acc = []\n",
    "dev_acc = []\n",
    "\n",
    "# Calcuate the number of parameter updates\n",
    "step = 1\n",
    "\n",
    "# Iterative training\n",
    "for epoch in range(max_iter):\n",
    "    # Random shuffle at the begging of each epoch\n",
    "    X_train, Y_train = _shuffle(X_train, Y_train)\n",
    "        \n",
    "    # Mini-batch training\n",
    "    for idx in range(int(np.floor(train_size / batch_size))):\n",
    "        X = X_train[idx*batch_size:(idx+1)*batch_size]\n",
    "        Y = Y_train[idx*batch_size:(idx+1)*batch_size]\n",
    "\n",
    "        # Compute the gradient\n",
    "        w_grad, b_grad = _gradient(X, Y, w, b)\n",
    "            \n",
    "        # gradient descent update\n",
    "        # learning rate decay with time\n",
    "        w = w - learning_rate/np.sqrt(step) * w_grad\n",
    "        b = b - learning_rate/np.sqrt(step) * b_grad\n",
    "\n",
    "        step = step + 1\n",
    "            \n",
    "    # Compute loss and accuracy of training set and development set\n",
    "    y_train_pred = _f(X_train, w, b)\n",
    "    Y_train_pred = np.round(y_train_pred)\n",
    "    train_acc.append(_accuracy(Y_train_pred, Y_train))\n",
    "    train_loss.append(_cross_entropy_loss(y_train_pred, Y_train) / train_size)\n",
    "    \n",
    "    y_dev_pred = _f(X_dev, w, b)\n",
    "    Y_dev_pred = np.round(y_dev_pred)\n",
    "    dev_acc.append(_accuracy(Y_dev_pred, Y_dev))\n",
    "    dev_loss.append(_cross_entropy_loss(y_dev_pred, Y_dev) / dev_size)\n",
    "\n",
    "print('Training loss: {}'.format(train_loss[-1]))\n",
    "print('Development loss: {}'.format(dev_loss[-1]))\n",
    "print('Training accuracy: {}'.format(train_acc[-1]))\n",
    "print('Development accuracy: {}'.format(dev_acc[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RJuoQ_R2jUmX"
   },
   "source": [
    "### Plotting Loss and accuracy curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 545
    },
    "colab_type": "code",
    "id": "DH3AJtvHjVJ7",
    "outputId": "4048b9f9-b739-4497-d2fe-93950bd25160"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9bnH8c+Tyb6whbCGkLBvCkJAVKCsCqig0gq4tdWK1Gptb6tVW63WerW91Var1iLordUL7gsVRQVBKaCsUpB9NawhbFnI/tw/ziRMwoRMIMmZTJ736zUvMmeZeWZe+j2/+Z1zfj9RVYwxxoSuMLcLMMYYU7cs6I0xJsRZ0BtjTIizoDfGmBBnQW+MMSHOgt4YY0KcBb0xxoQ4C3rTqInILhEZ7XYdxtQlC3pjjAlxFvTG+CEit4rINhE5IiLvi0g773IRkT+LyCEROS4i60Skj3fdeBH5RkSyRWSviPzS3U9hjMOC3phKRGQk8BhwLdAW2A3M8a6+FBgGdAOaAZOBLO+6WcBtqpoA9AEW1mPZxlQp3O0CjAlC1wMvqupqABG5DzgqIqlAEZAA9AC+UtWNPvsVAb1E5GtVPQocrdeqjamCteiNOV07nFY8AKqag9Nqb6+qC4FngGeBgyIyQ0SaeDedBIwHdovIYhG5qJ7rNsYvC3pjTrcP6Fj2RETigERgL4CqPq2qA4DeOF04d3uXr1DViUAr4F3g9Xqu2xi/LOiNgQgRiS574AT0D0Wkn4hEAf8NfKmqu0RkoIhcKCIRQC6QD5SISKSIXC8iTVW1CDgBlLj2iYzxYUFvDMwDTvo8hgIPAG8B+4HOwBTvtk2AF3D633fjdOn8ybvuRmCXiJwApgM31FP9xpyR2MQjxhgT2qxFb4wxIc6C3hhjQpwFvTHGhDgLemOMCXFBeWdsy5YtNTU11e0yjDGmwVi1atVhVU3yty4ogz41NZWVK1e6XYYxxjQYIrK7qnXWdWOMMSHOgt4YY0KcBb0xxoS4oOyjN8aYmioqKiIjI4P8/Hy3S6lT0dHRJCcnExEREfA+FvTGmJCQkZFBQkICqampiIjb5dQJVSUrK4uMjAzS0tIC3s+6bowxISE/P5/ExMSQDXkAESExMbHGv1os6I0xISOUQ77M2XzG0An64gL491Owe5nblRhjTFAJnaAvLYHlz8P8+6G01O1qjDGNzLFjx3juuedqvN/48eM5duxYHVR0SugEfWQsjHoQ9q2G9W+5XY0xppGpKuhLSs480di8efNo1qxZXZUFBBj0IjJWRDaLyDYRudfP+okisk5E1orIShEZUmm9R0TWiMi/aqtwv86fDG37woKHoehknb6VMcb4uvfee9m+fTv9+vVj4MCBjBgxguuuu47zzjsPgKuuuooBAwbQu3dvZsyYUb5famoqhw8fZteuXfTs2ZNbb72V3r17c+mll3LyZO3kWLWXV4qIB2fG+zFABrBCRN5X1W98NlsAvK+qKiLn48y52cNn/V3ARpxp2OpOWBhc+nv4x5Ww/G8w9L/q9O2MMcHp4bkb+GbfiVp9zV7tmvDbK3tXuf7xxx9n/fr1rF27lkWLFnH55Zezfv368ssgX3zxRVq0aMHJkycZOHAgkyZNIjExscJrbN26ldmzZ/PCCy9w7bXX8tZbb3HDDec+I2UgLfpBwDZV3aGqhcAcYKLvBqqao6fmJIwDyucnFJFk4HJg5jlXG4i0YdBtHHzxJORk1stbGmNMZYMGDapwrfvTTz9N3759GTx4MN9++y1bt249bZ+0tDT69esHwIABA9i1a1et1BLIDVPtgW99nmcAF1beSESuBh4DWuEEe5m/APcACWd6ExGZBkwDSElJCaCsMxjzO3huMCx+HC5/4txeyxjT4Jyp5V1f4uLiyv9etGgRn376KcuWLSM2Npbhw4f7vRY+Kiqq/G+Px1NrXTeBtOj9XbR52oziqvqOqvYArgIeARCRK4BDqrqqujdR1Rmqmq6q6UlJfodUDlxSN0i/GVa+BJmbz+21jDEmAAkJCWRnZ/tdd/z4cZo3b05sbCybNm1i+fLl9VpbIEGfAXTweZ4M7KtqY1X9HOgsIi2BS4AJIrILp8tnpIi8cvbl1sDweyEyDj55sF7ezhjTuCUmJnLJJZfQp08f7r777grrxo4dS3FxMeeffz4PPPAAgwcPrtfa5FTXehUbiIQDW4BRwF5gBXCdqm7w2aYLsN17MrY/MBdI9um3R0SGA79U1SuqKyo9PV1rZeKRJX+BT38LN70Pnb5z7q9njAlaGzdupGfPnm6XUS/8fVYRWaWq6f62r7ZFr6rFwB3AfJwrZ15X1Q0iMl1Epns3mwSsF5G1OFfoTNbqjiD14cLp0DQFPv61c0OVMcY0QgGNXqmq84B5lZY97/P3H4A/VPMai4BFNa7wXEREw+jfwlu3wLrXoN919fr2xhgTDELnztiq9JkE7QfAgkegMM/taowxpt6FftCLwGX/Ddn7YNkzbldjjDH1LvSDHiBlMPSc4JyczT7gdjXGGFOvGkfQA4x+CEoK4bNH3a7EGGPqVeMJ+sTOMOhWWPMKHNxQ/fbGGHOOHnroIf70pz+5XUYjCnqAYXdDVBP4+AG3KzHGmHrTuII+tgV85x7YvgC2fep2NcaYEPToo4/SvXt3Ro8ezebNzhAs27dvZ+zYsQwYMIChQ4eyadMmjh8/TmpqKqXeiZLy8vLo0KEDRUVFtV5TQNfRh5SBt8JXLzit+k4jIMzjdkXGmNr24b1w4D+1+5ptzoNxj59xk1WrVjFnzhzWrFlDcXEx/fv3Z8CAAUybNo3nn3+erl278uWXX3L77bezcOFC+vbty+LFixkxYgRz587lsssuIyIionbrpjEGfXgkjHkYXr8J1vwTBvzA7YqMMSHiiy++4OqrryY2NhaACRMmkJ+fz9KlS/ne975Xvl1BQQEAkydP5rXXXmPEiBHMmTOH22+/vU7qanxBD86llh0Gw8JHnRuqos44grIxpqGppuVdl0QqDvhbWlpKs2bNWLt27WnbTpgwgfvuu48jR46watUqRo4cWSc1Na4++jIicNmjkHsI/v2029UYY0LEsGHDeOeddzh58iTZ2dnMnTuX2NhY0tLSeOONNwBQVb7++msA4uPjGTRoEHfddRdXXHEFHk/ddCU3zqAHSE53WvNL/wrH97pdjTEmBPTv35/JkyfTr18/Jk2axNChQwF49dVXmTVrFn379qV3796899575ftMnjyZV155hcmTJ9dZXdUOU+yGWhumuDpHd8MzA53Av/pvdf9+xpg6Y8MUn8MwxSGteUcYPB2+ng37v3a7GmOMqRONO+gBhv4CYprD/F9DEP66McaYc2VBH90Uht8Hu76ALR+5XY0x5hwEY1d0bTubz2hBD5D+Q0js6txEVVL7d6UZY+pedHQ0WVlZIR32qkpWVhbR0dE12q9xXkdfmScCxvwO5kyFVf/rDH5mjGlQkpOTycjIIDMz0+1S6lR0dDTJyck12seCvkz3cZA6FBY9Budf63TpGGMajIiICNLS0twuIygF1HUjImNFZLOIbBORe/2snygi60RkrYisFJEh3uXRIvKViHwtIhtE5OHa/gC1RgQu/T3kHYEvnnS7GmOMqTXVBr2IeIBngXFAL2CqiPSqtNkCoK+q9gNuBmZ6lxcAI1W1L9APGCsig2ur+FrXrh/0nQLL/+ZcY2+MMSEgkBb9IGCbqu5Q1UJgDjDRdwNVzdFTZ0DiAPUuV1XN8S6P8D6C+0zJyN84rfuFj7hdiTHG1IpAgr498K3P8wzvsgpE5GoR2QR8gNOqL1vuEZG1wCHgE1X90t+biMg0b7fPSldPpjRNhovugP+8ARmr3KvDGGNqSSBBL36WndYqV9V3VLUHcBXwiM/yEm+XTjIwSET6+HsTVZ2hqumqmp6UlBRY9XVlyM8grhV8bDdRGWMavkCCPgPo4PM8GdhX1caq+jnQWURaVlp+DFgEjK15mfUsKgFG3A97lsHGuW5XY4wx5ySQoF8BdBWRNBGJBKYA7/tuICJdxDsIs4j0ByKBLBFJEpFm3uUxwGhgU21+gDpzwY2Q1BM+/S0UF7pdjTHGnLVqg15Vi4E7gPnARuB1Vd0gItNFZLp3s0nAem9f/LPAZO/J2bbAZyKyDueA8Ymq/qsuPkit84TDpY/AkR2wcpbb1RhjzFlr3MMUV0cV/nk17F8LP13jDH5mjDFByIYpPltlN1GdPAaf/8ntaowx5qxY0FenTR+44Ab48u9ON44xxjQwFvSBGPkb8ETCpw+5XYkxxtSYBX0gEtrAJXfBN+/BnuVuV2OMMTViQR+oi++AhLY2E5UxpsGxoA9UZJzThbN3JWx42+1qjDEmYBb0NdF3KrQ+z+mrL8p3uxpjjAmIBX1NhHngst/DsT3w1d/drsYYYwJiQV9TnYZD18vg8ycgN8vtaowxploW9GdjzO+gMAcW/8HtSowxploW9GejVQ8Y8H1nDJzD29yuxhhjzsiC/mwNvx/CY+CTB92uxBhjzsiC/mzFJ8HQn8PmD2DXErerMcaYKlnQn4vBt0PTDjD/figtdbsaY4zxy4L+XETEwKgHYf/XzhyzxhgThCzoz1Wf70K7C2DB76DopNvVGGPMaSzoz1VYGFz6KJzIgGXPul2NMcacxoK+NqReAj2ugCV/hpxDbldjjDEVWNDXltEPQ3E+fPbfbldijDEVBBT0IjJWRDaLyDYRudfP+okisk5E1orIShEZ4l3eQUQ+E5GNIrJBRO6q7Q8QNFp2gYE/gtX/gEMb3a7GGGPKVRv0IuIBngXGAb2AqSLSq9JmC4C+qtoPuBmY6V1eDPxCVXsCg4Gf+Nk3dHznVxCZYDdRGWOCSiAt+kHANlXdoaqFwBxgou8GqpqjWj4bRxyg3uX7VXW19+9sYCPQvraKDzqxLWDYL2Hrx7D6ZberMcYYILCgbw986/M8Az9hLSJXi8gm4AOcVn3l9anABcCX/t5ERKZ5u31WZmZmBlBWkLrwNmeEy/fvhIW/txupjDGuCyToxc+y0+bSU9V3VLUHcBXwSIUXEIkH3gJ+pqon/L2Jqs5Q1XRVTU9KSgqgrCAVHgXXvwn9b4LP/wfeusWurzfGuCqQoM8AOvg8Twb2VbWxqn4OdBaRlgAiEoET8q+qauOYg88TAVc+7QxnvOEd+MeVkNOAf6UYYxq0QIJ+BdBVRNJEJBKYArzvu4GIdBER8f7dH4gEsrzLZgEbVfXJ2i09yInAJXfB5H/CgfUwcyQc2uR2VcaYRqjaoFfVYuAOYD7OydTXVXWDiEwXkenezSYB60VkLc4VOpO9J2cvAW4ERnovvVwrIuPr5JMEq55Xwg/nQXEBzBoD2xe6XZExppGRUxfLBI/09HRduXKl22XUruMZ8H+TnWvsL38C0n/odkXGmBAiIqtUNd3fOrsztr40TYabP4Iuo+BfP4P5v4bSErerMsY0Ahb09SkqAabMhkG3wbJn4LUboTDX7aqMMSHOgr6+ecJh/B9h3B9hy4fw0jg4sd/tqowxIcyC3i0X3gZTX4Os7fDCSNi/zu2KjDEhyoLeTd0uhZvng4TBi2Nhy3y3KzLGhCALere16QO3LoCWXWH2FFj+PAThlVDGmIbLgj4YJLRxrrXvPh4++hXMuxtKit2uyhgTIizog0VkHFz7T+du2hUvOK37fL/DAhljTI1Y0AeTsDBnfJwrn4Idnzn99se+rX4/Y4w5Awv6YDTgB84ImMcznCty9q5yuyJjTANmQR+sOo+AH30CETHw0uXwzXtuV2SMaaAs6INZUne4dSG0PR9evwmW/NmuyDHG1JgFfbCLawk3vQ99JsGnDzkzVxUXul2VMaYBCXe7ABOAiGiYNAsSu8DiP8Cx3XDtyxDT3O3KjDENgLXoGwoRGHE/XD0D9iyHmWPgyA63qzLGNAAW9A1N38lw03uQlwUvjILdy9yuyBgT5CzoG6KOF8OPPoXYFvDyBFj3htsVGWOCmAV9Q5XYGW75BDpcCG//CBY9blfkGGP8sqBvyGJbwA1vQ7/rYdFj8PatUJTvdlXGmCATUNCLyFgR2Swi20TkXj/rJ4rIOu/k3ytFZIjPuhdF5JCIrK/Nwo1XeCRMfBZGPQj/eQNengi5h92uyhgTRKoNehHxAM8C44BewFQR6VVpswVAX1XtB9wMzPRZ97/A2FqpthpHcwvJLWiEoz6KwNBfwPf+F/avhZmjIHOL21UZY4JEIC36QcA2Vd2hqoXAHGCi7waqmqNa3kEcB6jPus+BI7VUb5WO5xUx+snFPLVga12/VfDqfTX84ANnHtpZo2HHYrcrMsYEgUCCvj3gO4RihndZBSJytYhsAj7AadXXiIhM83b7rMzMzKzp7jSNjWB0z9a8uGQnWw5m13j/kJGcDj9aAAnt4JVrYNmzdpLWmEYukKAXP8tOSw5VfUdVewBXAY/UtBBVnaGq6aqanpSUVNPdAbhnbHfiosL5zbvr0cYcbs07wi3zodtYmH8/vHYDnDzqdlXGGJcEEvQZQAef58nAvqo29nbVdBaRludYW40lxkdxz9jufLXzCO+trbLExiG6KUx+BS57DLZ8BH8fBntXu12VMcYFgQT9CqCriKSJSCQwBXjfdwMR6SIi4v27PxAJZNV2sYGYMjCFvslN+f0HGzmRX+RGCcFDBC663ZmAXBVmXQpf/t26coxpZKoNelUtBu4A5gMbgddVdYOITBeR6d7NJgHrRWQtzhU6k8tOzorIbGAZ0F1EMkTklrr4IGU8YcIjV/UhK7eAJz+2K08Ap9/+ts+hyyj48B544/uQf9ztqowx9USCsS87PT1dV65ceU6v8Zt3/8P/fbmHuXcOoXe7prVUWQNXWgrL/gqfPgzNUuDaf0Dbvm5XZYypBSKySlXT/a0L2Ttj7760B81iI3ng3fWUlgbfwcwVYWHO5OM/nAfFBc4ImCtmWVeOMSEuZIO+aWwE947rweo9x3hzdYbb5QSXlMEwfQmkDYUP/gve+hEUNOJLUo0JcSEb9ADf7Z/MgI7NefzDTRzLs1mZKohLhOvecIZO2PA2zBgOB2yUCmNCUUgHfViY8MjEPhzLK+R/5m92u5zgExbmDJ3w/blQkOMMnbD6ZevKMSbEhHTQA/Rq14TvX5zK/321h6+/PeZ2OcEpdQhM/8Lp0nn/TnhnujOMgjEmJIR80AP8fEw3WsZH8cB76ymxE7P+xbdyhjwefj+sew1eGAmHNrldlTGmFjSKoG8SHcFvLu/JuozjzP5qj9vlBK8wDwz/Fdz0rneqwhGwdrbbVRljzlGjCHqACX3bMbhTC/5n/maycgrcLie4dRruXJXTfgC8Ox3e+wkU5rldlTHmLDWaoBdxTszmFhTzh4+sS6JaCW3gxndh2N2w5lXnRO3hRjwEtDENWKMJeoCurRO4ZUgar6/MYNXuOh8iv+HzhMPI38ANb0LOQfj7d2wicmMaoEYV9AA/HdWVtk2j+c27GyguKXW7nIahy2i47Qtoe74zEfncn9nctMY0II0u6OOiwnngil5s3H+Cfy7f7XY5DUfT9s719pf8DFa95MxglbXd7aqMMQFodEEPMK5PG4Z2bcmTH2/h0AlrmQbMEwFjHobrXofjGU5XzoZ33K7KGFONRhn0IsLDE3pTUFzKYx/aidka63aZ05XTqge88QOYd7czSJoxJig1yqAH6JQUz7RhnXhnzV6W73BljpSGrVkH+ME8uOgO+GoGvHgZHN3ldlXGGD8abdAD/GREF9o3i+HB99ZTZCdmay48Ei57FCa/Clk74PlhsPFfbldljKmkUQd9TKSHhyb0ZsvBHF769063y2m4el4B0z+HxE7w2vXw0f1QbKOFGhMsGnXQA4zp1ZpRPVrxl0+3sv/4SbfLabiapzpz0w66DZY/Cy+Ng2M23IQxwaDRBz3AQxN6U1Kq/P6DjW6X0rCFR8H4P8L3/gGHt8DzQ2HzR25XZUyjFx7IRiIyFngK8AAzVfXxSusnAo8ApUAx8DNVXRLIvsGgQ4tYbh/ehT9/uoWpAw8zpGtLt0tq2HpfBW3OcyYhnz0Z+l0PLdIgIg4iKz3Kl8VCZLzztycSRNz+FMaEjGonBxcRD7AFGANkACuAqar6jc828UCuqqqInA+8rqo9AtnXn9qYHLym8otKuOwvn+MR4cOfDSUq3FOv7x+SivLh4187k5mU1KDPXjze0I/1Hgx8DgJlB4QI77oK2/k5aJTvGwvh0XYAMSHrTJODB9KiHwRsU9Ud3hebA0wEysNaVXN8to8DNNB9g0V0hIeHJ/TmBy+tYOYXO/nJiC5ul9TwRUTD5U84j+JCKMp1JjQpzIPCHCjK8z7P8S7LPcM2uc7Qycf2eJd59ympwfX7ca2g26XQbZwzQmdUfF19cmOCSiBB3x741ud5BnBh5Y1E5GrgMaAVcHlN9vXuPw2YBpCSkhJAWbVvePdWjO3dhr8u3MqEvu3o0CLWlTpCUnik84hpXruvW1LkPUD4HBDKn/scQApz4OB6+GYurHnF6R5KHQrdxzk3gDVz5785Y+pDIEHv77fuaf09qvoO8I6IDMPprx8d6L7e/WcAM8DpugmgrjrxwJW9WPxEJr/71ze8cJPfX0EmmHgiIKaZ8whESRHsWQ5bPnIe837pPFr1gm5jnUdyujMJizEhIpCgzwA6+DxPBvZVtbGqfi4inUWkZU33DQbtm8Vw56gu/PGjzSzcdJCRPVq7XZKpTZ4ISBvqPC57FA5vg63zYfOHsPRpWPIkxLSArpdC97HQeSREN3W7amPOSSAnY8NxTqiOAvbinFC9TlU3+GzTBdjuPRnbH5iLE+qe6vb1x42Tsb4Ki0sZ99TnFJUoH/98GNER1rprFE4eg+0LnZb+1o/h5FEIC4eOFzv9+t0ug8TObldpjF9nOhlbbdB7X2A88Bec4H5RVR8VkekAqvq8iPwKuAkoAk4Cd/tcXnnavtW9n9tBD7B022Gum/kld43qys/HdHO1FuOC0hLIWOG09LfMh0zvPRaJXZ2Wfrex0GGwMzmLMUHgnIO+vgVD0APcOXsN8zcc4JOfD6NjYpzb5Rg3Hd3lBP6Wj2DXEudy0eim0GWME/pdRkFsC7erNI2YBf1ZOngin1FPLCY9tTkv/WAgYtdgG4CCbNj+mRP8W+dDbqZz7X/KYKd7p9tYaNnNrtk39cqC/hzM/GIHv/9gI8/fMICxfdq4XY4JNqWlsG/1qat4DvzHWd489VS/fsdLnEtLjalDFvTnoKiklCueXkJOQTGf/NcwYiOtT9acwfGMU108OxY7N3RFJkCXkd4unjEQn+R2lSYEWdCfo692HuHavy/j9uGduWdsD7fLMQ1FYS7s/Nzb2p8P2fsBgfYDoOsYJ/TbXQBhNragOXcW9LXgv15fy9yv9/HhXcPo0spunTc1pAr7v/ZeuvkJ7F0FKMQmQpfRTuh3HglxiW5XahooC/pakJldwMgnFnF+clNeueVCOzFrzk3uYeea/a2fwPYFzjg+vq39rmOgrbX2TeAs6GvJy8t28eB7G/jr1Au4sm87t8sxoaK0BPathW2fODdq7V2N09pv6Vy2aa19EwAL+lpSUqpMeGYJh3MKWPCL4cRH2YlZUwfO2Nq/FLqOtta+OY0FfS1aveco1zy3lFuHpvHry3u5XY4JdaUlsG+NE/rbPvHf2rebtQznPh698dE/pTlTBnbgxX/v4rsDOtC9TYLbJZlQFuZxRtNMTocR9/m09j92wn/da4A467uMsda+8cta9GfhSG4hI59YRLfWCbw2bbCdmDXuqK613/VSp2/fjdZ+aQnkH3cGhss74vx78sjpfxflQevekHIxdBgE0U3qv9YQYV03dWD2V3u47+3/8OS1fbmmf7Lb5RjjtPa3LXBCf9sCJ0wlzOnbP9vWvqozaYvfsD5DiJ88RhVTTwDiTEAT09yZUD5zM2iJU2vrPpByEXS8yAn/BBsmPFAW9HWgtFS55m9LyTiax4JfDKdpTITbJRlzSnWt/c6jICLGTyv76OnBXVpU9ftEJkBsc2cM/5jmzq+HKv/2hnt0s4oHm4IcZ6TQPctg91LIWAnFJ511LTo5gd/xIucA0KKTjSFUBQv6OrJ+73GufGYJ378olYcm9Ha7HGOq5q+17ys82n8oVxfcnjpo4JQUOTeX7V7qhP+eZc4BByC+tTN4XFn4t+5js4F5WdDXoQfeXc+rX+7m/TuG0Ke9zURkGoDSEji4wWkZlwV3ZBDPj1xaCoc3e1v83uA/7p2KOqqJ07dfFv7tBziT0jdCFvR16HheESOfWERKYixvTb+YsDD7WWlMnTv27amunj3LT00M44mEdv1PdfV0uDDw+YQbOAv6Ovbmqgx++cbX/GHSeUwemOJ2OcY0PnlHnMDfs9Rp9e9fC6XFgHiv6vE5wdukrdvV1gkL+jqmqlz792VsO5TDwl8Mp3mcjT1ujKsKc52B43Yvc8L/2xVQlOusa9bRmQc45SLn38QuIXGC14K+Hmw6cILLn17CtekdeOya89wuxxjjq6QIDqw71ce/Z5l3aAkgLgmSBzpX9DRPheZpzr/NUhrUhDHnfGesiIwFnsKZ4Humqj5eaf31wK+8T3OAH6vq1951dwG3AgK8oKp/OatPEeR6tGnCDy5O5cV/72TywA7069A4+gWNaRA8Ec6J2vYD4OI7nPsDDm891dWzb41zNVJJgc9OAk2TveHfseJBoHmqM8R0A/klUG2LXkQ8wBZgDJABrACmquo3PttcDGxU1aMiMg54SFUvFJE+wBxgEFAIfIRzENh6pvdsiC16gOz8IkY9sZjWTaJ568cXExlut6Eb02CUlkLOQWcieH+PnAMVt4+MPxX6lR/NUpybwerRubboBwHbVHWH98XmABOB8qBX1aU+2y8Hym4V7QksV9U8776LgauBP9b0QzQECdER/Pryntw1Zy09H/yI5OYxdEyMIy0xltSWcaQmxpHaMo7k5jFEeOwgYExQCQtzTtQ2aeucuK2sMA+O7Tn9AJC13fk1UHaTFwACTdpVfSCIS6rXXwOBBH174Fuf5xnAhWfY/hbgQ+/f64FHRSQROAmMBxpeU70GJvRtR4fBQqoAABD3SURBVGxkOOsyjrHzcC67snJZvfsoOQXF5dt4woQOZQeBlnGkJsbSsWUcaYnOQSDcDgLGBJ/IWGjVw3lUpgo5h/z/Etj+GWTvq7h9RGzVB4Gk7rVeeiBB7++w47e/R0RG4AT9EABV3SgifwA+wem7/xoormLfacA0gJSUhnuJoogwpldrxvQ6NUaHqpKVW8iuw7nsPJzL7qw8dmblsutwLit3HSG3sKR82/AwIbl5TPkvgLSWcXRMjCWtZRztm9lBwJigJOKMy5PQGlL8tIOLTjrX/vs7EOxYfOqKoNhEuGdHrZcXSNBnAB18nicD+ypvJCLnAzOBcaqaVbZcVWcBs7zb/Lf39U6jqjOAGeD00QdYf4MgIrSMj6JlfBTpqRVHElRVDucUsssb/M6/eew8nMuKnacfBDq0iCW1UldQWmIc7ZpF20HAmGAVEQNJ3ZxHZarOEBVHd0H+sTp5+0CCfgXQVUTSgL3AFOA63w1EJAV4G7hRVbdUWtdKVQ95t7kG8NP51XiJCEkJUSQlRDHQz0EgM6eAXYfzTjsQfLnzCHk+B4EIj9Chue8BIJburRMYlNbChlE2JpiJQHyS86gj1Qa9qhaLyB3AfJzLK19U1Q0iMt27/nngQSAReM4bKsU+Z3/f8vbRFwE/UdWjdfA5QpKI0CohmlYJ0QxK83MQyC5gV1ae0yWUlcvurFx2Hs5j+Y6s8oPA5ee15bFJ59Ek2kbXNKaxshumQlDZQeCt1Xv508ebadcsmmem9qevXdtvTMg60+WV1qkbgkSEVk2i+fHwzrx+20WUlsJ3n1/KzC92EIwHdmNM3bKgD3EDOjZn3k+HMrJHK37/wUZu+cdKjuQWul2WMaYeWdA3Ak1jI3j+hgE8PKE3S7YeZvxTX/DljqzqdzTGhAQL+kZCRPj+xam8ffvFxER6mPrCcp5esJWSUuvKMSbUWdA3Mn3aN2XunUOY0LcdT36yhRtnfcmhE/lul2WMqUMW9I1QfFQ4f57cjz9+93zW7DnGuKe+YPGWTLfLMsbUEQv6RkpEuDa9A+/fcQkt46P4/otf8YePNlFUUup2acaYWmZB38h1bZ3Ae3dcwtRBKfxt0XYm/30ZGUfz3C7LGFOLLOgN0REeHrvmPP469QK2HMxh/FNfMH/Dgep3NMY0CBb0ptyVfdvxwU+H0DExjtv+uYqH3t9AQXFJ9TsaY4KaBb2poGNiHG/++CJuviSN/126i2ueW8rOw7lul2WMOQcW9OY0UeEeHryyFzNvSmfvsZNc8fQXvLd2r9tlGWPOkgW9qdLoXq2Z99Oh9GzbhLvmrOWeN78mr9DvvDHGmCBmQW/OqF2zGOZMG8wdI7rwxqoMJj7zbzYfyHa7LGNMDVjQm2qFe8L45WXd+efNF3I0r4gJzyxh9ld7bCRMYxoIC3oTsCFdWzLvriEMTG3BfW//hztnryE7v8jtsowx1bCgNzXSKiGal28exN2XdefD9Qe4/OklrMuom3kujTG1w4Le1FhYmPCTEV14bdpgiktKmfS3pcxastO6cowJUhb05qylp7Zg3l1D+U63Vjzyr2+49eVVHLVJTYwJOgEFvYiMFZHNIrJNRO71s/56EVnnfSwVkb4+634uIhtEZL2IzBaR6Nr8AMZdzWIjeeGmATx4RS8WbznE+Ke/YMWuI26XZYzxUW3Qi4gHeBYYB/QCpopIr0qb7QS+o6rnA48AM7z7tgd+CqSrah/AA0ypvfJNMBARbh6Sxls/vpjI8DCmzFjOs59to9QmNTEmKATSoh8EbFPVHapaCMwBJvpuoKpLVfWo9+lyINlndTgQIyLhQCyw79zLNsHo/ORm/OvOIYw/ry3/M38zN734FYeybVITY9wWSNC3B771eZ7hXVaVW4APAVR1L/AnYA+wHziuqh+fXammIUiIjuDpKf14/JrzWLHrCOOfWsJH6/fbhOTGuCg8gG3EzzK/v8lFZARO0A/xPm+O0/pPA44Bb4jIDar6ip99pwHTAFJSUgIq3gQnEWHKoBQuSGnOHf+3mumvrAageWwEnZLi6ZwUR+ek+PK/U1rEEu6x6wKMqSuBBH0G0MHneTJ+ul9E5HxgJjBOVbO8i0cDO1U107vN28DFwGlBr6oz8Pbtp6enW+duCOjeJoG5dw5h2Y4sth/KYXtmLjsyc1i4KZPXV2aUbxfhEVJaxNI5KZ7OreLp1DKOzq3i6dwynqaxES5+AmNCQyBBvwLoKiJpwF6ck6nX+W4gIinA28CNqrrFZ9UeYLCIxAIngVHAytoo3DQM0REeRnRvxYjurSosP36yiB2ZTvhvz8wp//uzzYcoKjl1nG8ZH+lt+Z/6JdA5KZ72zWPwhPn7sWmMqazaoFfVYhG5A5iPc9XMi6q6QUSme9c/DzwIJALPiQhAsaqmq+qXIvImsBooBtbgbbWbxq1pTAQXpDTngpTmFZYXl5Ty7dGT3l8AOezwHgg+Wr+fo3mnhluIDA8jLTGOzq3i6NQyvvzfTklxJETbrwBjfEkw3s2Ynp6uK1daw99UdCS30NvyP3UA2J6Zy54jeZT4XMrZuklUefj7ngto1zSGMPsVYEKUiKxS1XR/6wLpujEmKLSIi6RFXAvSU1tUWF5YXMqeI7lsO5TLjsM5bD/kHATeW7uP7PxT4+fHRnro3iaBHm2a0KttAj3aNqFHmwT7BWBCngW9afAiw8Po0iqBLq0SKixXVQ7nFHpb/jlsPZjDxv0n+GDdPmZ/deoAkNw8pkL492zbhJQWsXYOwIQMC3oTskSEpIQokhKiGNwpsXy5qrL/eD6bDpxg4/5sNu4/waYD2SzcdJCyHqCYCA/d2iQ44d/Gafn3aNuEpjHW+jcNj/XRG+OVX1RS3urfeOAEm/Zns/HACY75nARu3yyGnmXh3zaBnm2bkJoYZ61/4zrrozcmANERHs5Lbsp5yU3Ll6kqB08UsPHACafl7/0F8NnmzPITwNERYXRvXTH8e7ZpYvcAmKBhQW/MGYgIbZpG06ZpdIV7AfKLSth2KKe822fTgRN8svEgr608NVpIu6bR5Sd8e7ZtQs+2CaQmxtldwKbeWdAbcxaiIzz0ad+UPu0rtv4zswv4piz89zvnAD7fkkmxt/UfFR5G6ybRxEeFkxBd9oggPiqc+LLnUf6WRZAQ7TyPsAOFqSELemNqiYjQqkk0rZpEM9yn9V9QXML2Q7ne1v8JMrMLyM4vJrugmH3H8skuyCYnv5js/OLyA8KZRIWHkRAdUX6giI8K9x44Ki5LiI5wDhTeg0q870ElKtzOKzQiFvTG1LGocA+92jWhV7smZ9xOVSkoLuVEfhE5+cXkFDjh7zyKyCkodg4IBacv252TR05BsbNvQTGBXGPRpkk0KYmxdGwRS8fEWFIS4+jYIpbUxDg7vxBiLOiNCRIiQnSEh+gID5VuCagRVSW3sMR7sCjiRH5x+S+GnIIisvOLOXGyiIxjJ9mTlceiLZlkZhdUeI2mMRFO+HsPAh1bxDn/JsbRKiHK7jBuYCzojQkxIlLePQOBzdyZV1jMniN57M7KY3dWLruz8thzJI91Gcf5cP2BCkNMRIWHlR8AUlrEkdqy7IAQR/tmMUSG2zmEYGNBb4whNjLce2PY6d1LRSWl7Dt2kl1ZeezxHgR2H8ljT1YeS7YdJr+otHzbMIF2zWJOHQQSTx0QOibGEhdlkeMG+9aNMWcU4QmjY2IcHRPjgKQK68quNNrl/SVQ/qvgSN5pI46CM+x0Weu/Y2Is7Zs5w02XlKrzUD31t++yEuff0lKluIplpaoUl1R8jbJlperdz8+y0lIlwhNGl1bxdG+TUD4eUou4yHr8luuW3RlrjKkzx08WsScrj91HvN1BPn8fOJEf0EljXyLgEcETVukR6DLv8rAwIdxn2cnCErYeyqkw5WVSQpQz9EWbBLp7h8Ho0iqe6AhPLX9LtcPujDXGuKJpTMRpdxuXyS8q4eAJZ/L4MG8oh4edCuGwymHtDei6oqpk5hSwaX82mw9ks+lANpsPnuDlZbspKC711gmpLeOc8G/t3Ando00CHZrHBvUJagt6Y4wroiM83u6g4CAitEqIplVCNMO6neqiKilVdmXlOuHvvRluw74TfLj+QPkvkthID11bJ9CjdVnXj/NvYnyUS5+mIuu6McaYs5BXWMyWgzlsPlB2J3Q2mw9m++3+6d76VN9/19Z10/1jXTfGGFPLYiPD6dehGf06NCtfVtb9s/mAT/fPgWz+ubxS909iHD3aOt0/Zb8AUlrUXfePBb0xxtQS3+6foV2r6P45kM3mAyf4plL3T0yEhz7tm/D6bRfhnXu71ljQG2NMHfOECZ2T4umcFM/489qWL6/c/XOysKTWQx4CDHoRGQs8BXiAmar6eKX11wO/8j7NAX6sql+LSHfgNZ9NOwEPqupfzrlyY4xp4Px1/9SFaoNeRDzAs8AYIANYISLvq+o3PpvtBL6jqkdFZBwwA7hQVTcD/XxeZy/wTi1/BmOMMWcQyKAUg4BtqrpDVQuBOcBE3w1UdamqHvU+XQ4k+3mdUcB2Vd19LgUbY4ypmUCCvj3wrc/zDO+yqtwCfOhn+RRgdlU7icg0EVkpIiszMzMDKMsYY0wgAgl6f2cG/F58LyIjcIL+V5WWRwITgDeqehNVnaGq6aqanpSUVNVmxhhjaiiQk7EZQAef58nAvsobicj5wExgnKpmVVo9DlitqgfPtlBjjDFnJ5AW/Qqgq4ikeVvmU4D3fTcQkRTgbeBGVd3i5zWmcoZuG2OMMXWn2ha9qhaLyB3AfJzLK19U1Q0iMt27/nngQSAReM57DWhx2a24IhKLc8XObXXzEYwxxpyJjXVjjDEh4Exj3QRl0ItIJnC2l2G2BA7XYjkNmX0XFdn3UZF9H6eEwnfRUVX9XskSlEF/LkRkZVVHtcbGvouK7PuoyL6PU0L9u7BZfI0xJsRZ0BtjTIgLxaCf4XYBQcS+i4rs+6jIvo9TQvq7CLk+emOMMRWFYoveGGOMDwt6Y4wJcSET9CIyVkQ2i8g2EbnX7XrcJCIdROQzEdkoIhtE5C63a3KbiHhEZI2I/MvtWtwmIs1E5E0R2eT9b+Qit2tyk4j83Pv/yXoRmS0i0W7XVNtCIuh9JkcZB/QCpopIL3erclUx8AtV7QkMBn7SyL8PgLuAjW4XESSeAj5S1R5AXxrx9yIi7YGfAumq2gdnmJcp7lZV+0Ii6AlgcpTGRFX3q+pq79/ZOP8jn2kOgZAmIsnA5TijqzZqItIEGAbMAlDVQlU95m5VrgsHYkQkHIjFz+i8DV2oBH1NJ0dpNEQkFbgA+NLdSlz1F+AeoNTtQoJAJyATeMnblTVTROLcLsotqroX+BOwB9gPHFfVj92tqvaFStAHPDlKYyIi8cBbwM9U9YTb9bhBRK4ADqnqKrdrCRLhQH/gb6p6AZALNNpzWiLSHOfXfxrQDogTkRvcrar2hUrQBzQ5SmMiIhE4If+qqr7tdj0uugSYICK7cLr0RorIK+6W5KoMIENVy37hvYkT/I3VaGCnqmaqahHOvBoXu1xTrQuVoK92cpTGRJxJAWYBG1X1SbfrcZOq3qeqyaqaivPfxUJVDbkWW6BU9QDwrYh09y4aBXzjYklu2wMMFpFY7/83owjBk9OBTCUY9KqaHMXlstx0CXAj8B8RWetddr+qznOxJhM87gRe9TaKdgA/dLke16jqlyLyJrAa52q1NYTgcAg2BIIxxoS4UOm6McYYUwULemOMCXEW9MYYE+Is6I0xJsRZ0BtjTIizoDfGmBBnQW+MMSHu/wEvT/Jx52BOTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEICAYAAAC0+DhzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yV9fXA8c/JICEk7DDDCLJlExBlqIACgjgRcKMVd1Hrz1VbtdZW66haVOpArTJFUHAioigiSthbwkyYYYaRhIzz++N7wUsSSAhJnpvkvF+vvMh95nkucM/9blFVjDHGGH9BXgdgjDEm8FhyMMYYk4slB2OMMblYcjDGGJOLJQdjjDG5WHIwxhiTiyUHY4wxuVhyMOWaiHwvIvtEJMzrWIwJJJYcTLklIo2BnoACg0vwviEldS9jCsuSgynPbgTmA+8BNx3bKCIVReRFEdksIgdEZK6IVPTt6yEi80Rkv4gkisjNvu3fi8gf/K5xs4jM9XutInK3iKwD1vm2veK7RoqILBSRnn7HB4vIYyKyXkQO+vY3EJHXRORF/4cQkRkicl9xvEGm/LLkYMqzG4Fxvp9+IlLbt/0FoDNwHlAdeAjIFpGGwJfAf4BooAOw5DTudzlwDtDa93qB7xrVgfHARyIS7tv3ADAcuASoDNwCHAHeB4aLSBCAiNQE+gATTufBjcmPJQdTLolID6ARMFlVFwLrgWt9H7q3AKNUdauqZqnqPFVNB64DZqnqBFXNUNU9qno6yeGfqrpXVVMBVPVD3zUyVfVFIAxo4Tv2D8DjqrpWnaW+Y38FDuASAsAw4HtV3XmGb4kxJ7DkYMqrm4CZqrrb93q8b1tNIByXLHJqcJLtBZXo/0JE/iQiq31VV/uBKr7753ev94Hrfb9fD3xwBjEZkydrGDPljq/94BogWER2+DaHAVWBukAacBawNMepiUDXk1z2MBDh97pOHsccnwLZ177wMK4EsFJVs0VkHyB+9zoLWJHHdT4EVohIe6AV8MlJYjKm0KzkYMqjy4EsXN1/B99PK+BHXDvEWOAlEannaxg+19fVdRzQV0SuEZEQEakhIh1811wCXCkiESLSFLg1nxiigEwgGQgRkb/i2haOeRt4WkSaidNORGoAqGoSrr3iA+DjY9VUxhQlSw6mPLoJeFdVt6jqjmM/wGhcu8IjwHLcB/Be4DkgSFW34BqI/+TbvgRo77vmv4GjwE5ctc+4fGL4Gte4/RuwGVda8a92egmYDMwEUoB3gIp++98H2mJVSqaYiC32Y0zpIyK9cNVLjVU12+t4TNljJQdjShkRCQVGAW9bYjDFxZKDMaWIiLQC9uMazl/2OBxThlm1kjHGmFys5GCMMSaXMjHOoWbNmtq4cWOvwzDGmFJl4cKFu1U1Oq99ZSI5NG7cmPj4eK/DMMaYUkVENp9sn1UrGWOMycWSgzHGmFwsORhjjMmlTLQ55CUjI4OkpCTS0tK8DqXYhYeHExMTQ2hoqNehGGPKiDKbHJKSkoiKiqJx48aISP4nlFKqyp49e0hKSiI2NtbrcIwxZUSZrVZKS0ujRo0aZToxAIgINWrUKBclJGNMySmzyQEo84nhmPLynMaYklNmq5WMKQuOHM1k0oJEBKhTpSJ1qoRTp3I40VFhBAfZlwKvqCrZSpn+O7DkUIz279/P+PHjueuuu07rvEsuuYTx48dTtWrVYorMlAbzEnbz8NRlJO7NvZZPkECtqHBqVwmnTuUw6lR2v9etEk7tyi6B1KkSTkQF+y9elHYdTOPjhVuZHJ/Ixt2HiQoLoXLFUKLCQ6hSMZTKFUOpHB7q+z2EyuFuW5WKoVQOD/n994qhVKoQHNClfvuXU4z279/P66+/nis5ZGVlERwcfNLzvvjii+IOzQSwg2kZ/PPLNYz/ZQuxNSsxaWQ3mkRHsuNAGjtS3M9O3+87U9JYn3yYeQl7OJiemetaUeEhuRJG7cp+SaRKONUjKhBUhr8Bn6nMrGy+X5vMpPhEZq/ZRVa20jW2Ope2q8vB9ExSUjNJScvgQGoGiXuPcDAtkwOpGRzK4+/DX3CQHE8YORPK74nGtz+PY8JDT/4ZUhQsORSjRx55hPXr19OhQwdCQ0OJjIykbt26LFmyhFWrVnH55ZeTmJhIWloao0aNYuTIkcDv04EcOnSIAQMG0KNHD+bNm0f9+vX59NNPqVixYj53NqXV92t38djU5exISeO2nrE8cFELKlZwHwLRUWG0pcpJzz2cnnk8cWz3Sx47Drg/f9t5kOSD6WTnmIg5NFioFeVLGL5qK/+SSJ3K4dSqHEZYSPF+GAWazXsOMzk+kY/ik9h1MJ2akWHc1rMJ18TF0CQ6Mt/zM7OyOZTuEkVKqu/PtAxSUjP8fj9x+/YDqaSkZZKSmkF65qmX6qgQEkSViqFc0qYOT13Wpqge+7hykRyemrGSVdtSivSaretV5olLzz7lMc8++ywrVqxgyZIlfP/99wwcOJAVK1Yc73I6duxYqlevTmpqKl26dOGqq66iRo0aJ1xj3bp1TJgwgbfeeotrrrmGjz/+mOuvv75In8V478CRDJ7+fBVTFibRtFYkU+48j04Nq53WNSqFhXBWdCRnneKDKzMrm+RD6ccTxo4DaWz3K4ms2pbC7NW7SM3IynVudFQYPZrWpHfLWvRqHk2VimVvXE1aRhZfr9zBxF8T+XnDHoIELmhRi6FdGtC7ZS1CgwvehyckOIiqERWoGlGh0LHklUBSUjOOJ5ADqRk0rZV/oiqMcpEcAkXXrl1PGIvw6quvMm3aNAASExNZt25druQQGxtLhw5uDfvOnTuzadOmEovXlIxZq3by2LTl7Dl8lLsvPIt7ezcrtiqDkOAg6lapSN0qJy99qiopaZnHq7GOJY71yYf4fu0upi3eSnCQ0KVxNfq2qk3vlrUK9E06kK3alsLk+ESmLd7KgdQMGlSvyIMXN+fqzg2oUyXck5jCQ4MJDw2mVpQnty8fySG/b/glpVKlSsd///7775k1axY///wzERERXHDBBXmOVQgLCzv+e3BwMKmpuRsnTem07/BRnpyxkk+XbKNlnSjeuakLbWNOXm1UUkSEKr6G0xZ1TvxkyspWliTu49vVu5i9Zhd//3w1f/98NbE1K9G7ZS36tKxFXOPqVAgJ/F7yKWkZzFi6jUkLElmWdIAKwUH0b1OHYV0a0K1JjXLfDlMukoNXoqKiOHjwYJ77Dhw4QLVq1YiIiGDNmjXMnz+/hKMzXvpi+Xb++ukK9h/J4L6+zbjrgqal4gM1OEjo3Kg6nRtV56H+LUnad4Tv1uxi1updfPDzZt6Zu5GosBB6NY+md8taXNAimhqRYflfuISoKvGb9zHx10Q+X76NtIxsWtaJ4olLW3N5h/pUq1S4KqCyyJJDMapRowbdu3enTZs2VKxYkdq1ax/f179/f8aMGUO7du1o0aIF3bp18zBSU1KSD6bzxPQVfLF8B23qV+aDW8+hVd3KXodVaDHVIrjh3MbccG5jDqdn8lPCbmav2cW3a3bx+fLtiECnhtVcqaJVLVrUjvKk+2bywXSmLkpiUnwiG5IPExkWwhUdYxjWpQHtYqoEdJdSr5SJNaTj4uI052I/q1evplWrVh5FVPLK2/OWNqrK9KXbeHL6Sg6nZzGqbzNu79WEkNNo4CxNsrOVldtS+HbNTr5dvYvlWw8AUL9qRXq3rEXvVrU4t0mNYu2OmZWt/PBbMpMWJDJr9U4ys5W4RtUY2qUBA9vVtTEggIgsVNW4vPbZu2NMMduZksafp61g1uqddGxYleevbkdTr1oZS0hQkNA2pgptY6pwX9/m7ExJ4ztfiWLKwiQ+mL+ZiqHB9GhWkz4ta3Fhy1rUrlw0Db+Je4/wUXwiHy1MYvuBNGpUqsAtPWK5Ji6mzL/vRcmSgzHFRFWZsjCJpz9bRXpmNo8PbMWI7rFlesqFk6ldOZxhXRsyrGtD0jKymL9hj6t+Wr2Lb1btBKBt/SrHq5/a1KtyWg3C6ZlZzFy5k0kLEpmbsBsROL95NH8d1Jo+rWqXivacQGPJwZhisHV/Ko9NXc6c35Lp2rg6z13djtialfI/sRwIDw3mgha1uKBFLZ4arKzdefB476dXZ6/jlW/XUSsqzFU/taxFj2Y1T1oFtHbHQSYu2MK0xVvZfySD+lUrcn/f5lwdF0P9qjZY9ExYcjCmCKkq43/dwj+/WEO2Kk8NPpsbujUq990iT0ZEaFmnMi3rVObuC5uy51A6c35Ldg3ay7YzcUEiFUKCOLdJDfq0csmiakSF411QlyTuJzRYuPhs1wW1+1k17b0uIgVKDiLSH3gFCAbeVtVnc+yvAnwINPRd8wVVfde3737gD4ACy4ERqpomIs8DlwJHgfW+7ftFpDGwGljru/x8Vb3jTB7SmJKQuPcID3+8jHnr93DeWTV47qp2NKge4XVYpUqNyDCu7BTDlZ1iyMjKZsHGvXy7Zhffrt7JXz9dyV8/XUmF4CCOZmXTrFYkjw9sxZWdYqhuXVCLXL7JQUSCgdeAi4AkYIGITFfVVX6H3Q2sUtVLRSQaWCsi44Bo4I9Aa1VNFZHJwDDgPeAb4FFVzRSR54BHgYd911uvqh2K5hGNKV7Z2cr/ft7Ec1+tJThI+McVbRnetYF1jzxDocFBnNe0Juc1rclfBrVmQ/IhZq/Zxbb9aQxqX5eODarae1yMClJy6AokqOoGABGZCFwG+CcHBaLE/U1FAnuBY1MShgAVRSQDiAC2AajqTL/z5wNXn8FzlApPPvkkkZGRPPjgg16HYorIxt2HeWjKUhZs2sf5zaP555VtqWd13cWiSXRkqZ+mozQpSHKoDyT6vU4CzslxzGhgOu6DPwoYqqrZwFYReQHYAqQCM3MkhWNuASb5vY4VkcVACvC4qv6Y8wQRGQmMBGjYsGEBHsOYopOVrYydu5EXZq4lLCSI569ux9WdY+ybrCkzCtK/K69/7TlHzvUDlgD1gA7AaBGpLCLVcKWMWN++SiJywpSiIvJnXCljnG/TdqChqnYEHgDGi0iuIaSq+qaqxqlqXHR0dAEewxvPPPMMLVq0oG/fvqxd65pR1q9fT//+/encuTM9e/ZkzZo1HDhwgMaNG5Od7abpPXLkCA0aNCAjI8PL8E0e1u08yFVvzOOZL1bTs1k03zxwPkPirBrJlC0FKTkkAQ38XsfgqxryMwJ4Vt1w6wQR2Qi0BBoBG1U1GUBEpgLn4RqvEZGbgEFAH9+5qGo6kO77faGIrAeaAycOgT4dXz4CO5YX+vQ81WkLA5495SELFy5k4sSJLF68mMzMTDp16kTnzp0ZOXIkY8aMoVmzZvzyyy/cddddzJ49m/bt2zNnzhwuvPBCZsyYQb9+/QgNLXvTIpdWmVnZ/PeHDbwyax2VwoJ5ZVgHBrevZ0nBlEkFSQ4LgGYiEgtsxTUoX5vjmC1AH+BHEakNtAA24Eod3UQkAlet1Affh7yvB9TDwPmqeuTYhXwN2ntVNUtEmgDNfNcqdX788UeuuOIKIiJcj5XBgweTlpbGvHnzGDJkyPHj0tPTARg6dCiTJk3iwgsvZOLEiae9vGhZkZaRRfLBdCLDQogKDwmIKSZWb0/h/6YsZcXWFAa2rcuTg88mOipwJpQzpqjlmxx8vYnuAb7GdWUdq6orReQO3/4xwNPAeyKyHJcQHlbV3cBuEZkCLMJVHS0G3vRdejQQBnzj++Z1rMtqL+BvIpIJZAF3qOreM3rKfL7hF6ec3yqzs7OpWrUqS5YsyXXs4MGDefTRR9m7dy8LFy6kd+/eJRVmwFiWtJ/bP1jI9gO/T18eHhpEVHgoUb5kERUeejxxRPpeR+V4HRkWQmW/1xGhwYXq/340M5vXvkvgte8SqBoRyhvXdWJA27pF+cjGBCSbeK8YLVq0iJtvvplffvnleLXS7bffztSpU7n//vsZMmQIqsqyZcto3749AEOGDCE8PJyoqChef/31At8rEJ73TE1ZmMRj05YTHRnGXReeRXqGW2bxYJpbjzclLZNDab+/PnjsdT5r9QKI4BJKmC95hPuSie91VLjb559cRODf3/zGmh0HubxDPf566dnWn96UKTbxnkc6derE0KFD6dChA40aNaJnz54AjBs3jjvvvJO///3vZGRkMGzYsOPJYejQoQwZMoTvv//ew8hLVkZWNs98vpr35m3i3CY1GH1tx9NaAyA7Wzl81JcsfMnkYNqJrw+l+ZKLX7LZe/goW/Yc8W3PIC0j95q9taLCeOvGOC5qXTuPOxtTdlnJoYworc+7+1A6d49bxC8b93Jrj1geHdDSszaGo5nZHPaVSA6mZ3A4PYtWdaOICrdOAaZsspKDCUjLkvZzxwcL2XP4KP8e2p4rOsZ4Gk+FkCAqhFSw1cCMwZKD8cjHC5N41Ne+8PGd59GmvvdrJxtjflemk4Oqlos+6KWpavBM2xeMMSWjzCaH8PBw9uzZQ40aNcp0glBV9uzZQ3h40ayiVZwCqX3BGHNqZTY5xMTEkJSURHJystehFLvw8HBiYrytr89PoLUvGGNOrcwmh9DQUGJjY70Ow2DtC8aURmU2ORjvWfuCMaWXJQdTLKx9wZjSzZJDAPhy+Xae/WoNg9vX44ZujahVOfAbl09ledIBbv8g3toXjCnF7KucxzbuPsyDHy3lcHoWo79LoPtzs3lg0hJWbD3gdWiF8vHCJK4aMw8R4eM7z7PEYEwpZSUHD6VlZHHXuEWEhgQx/Z7upGdm8/68TUyOT2Tq4q10ja3OLd1juah1bYILMaNoSbL2BWPKljI7t1Jp8Pgny/lw/hbG3hxH75a/T+x2IDWDj+ITefenTWzdn0qD6hW5+bxYromLCch5fqx9wZjS6VRzK1ly8Mhny7Zxz/jFjOzVhMcuyXvCvMysbL5ZtZN35m4kfvM+IsNCuCauATef15iGNSJKOOK8+bcvPHtVW6tGMqYUseQQYDbtPsyg/8ylWe1IJt9+LqEF+Ja9NHE/7/60kc+WbSdblYta1+aW7rF0ja3u2Qhw//EL/72hs41fMKaUseQQQNIzs7jqjXkk7k3li1E9qV+14mmdv+NAGh/M38S4X7aw/0gGbepX5pbusQxqV48KISVTlePfvtCtSXVeu7aTtS8YUwpZcgggT3y6gvd/3nzGC8ikHs1i2uKtjP1pIwm7DhEdFcaN3Rpx7TkNi/WD2r994ZbusTx2ibUvGFNaWXIIEF8u386d4xbxhx6xPD6odZFcU1X5Yd1uxs7dyJzfkqkQEsQVHepzS49YWtSJKpJ7HOPfvvDPK9tyZSdrXzCmNLPFfgLAlj1HeGjKMto3qMpD/VsW2XVFhPObR3N+82gSdh1k7E+bmLooiUnxifRoWpNbe8RyfvNogs6wK6zNj2RM+VKg+gAR6S8ia0UkQUQeyWN/FRGZISJLRWSliIzw23e/b9sKEZkgIuG+7dVF5BsRWef7s5rfOY/67rVWRPoVxYN6KT0zi3smLEIERg/vWGxtA01rRfGPK9ry8yN9eKh/CxJ2HWLEewvo+9IcPvh5E0eOZp72NTOysnlqxkr+9NFSOjWsyvR7ultiMKYcyLdaSUSCgd+Ai4AkYAEwXFVX+R3zGFBFVR8WkWhgLVAHiAbmAq1VNVVEJgNfqOp7IvIvYK+qPutLONV857cGJgBdgXrALKC5qmadLMZAr1Z6asZK3v1pE/+9oTP9zq5TYvfNyMrmi+XbGTt3I0uTDlA5PITh5zTkpnMbU68ADeF7DqVz9/hFzN9g7QvGlEVnWq3UFUhQ1Q2+i00ELgNW+R2jQJS4PpWRwF7g2NfUEKCiiGQAEcA23/bLgAt8v78PfA887Ns+UVXTgY0ikuCL4ecCxBpwvlqxg3d/2sSI7o1LNDEAhAYHcVmH+gxuX49FW/bxztyNvPXDBt7+cSMD2tThlh6xdGpYLc9z/dsXXrqmvbUvGFPOFCQ51AcS/V4nAefkOGY0MB33wR8FDFXVbGCriLwAbAFSgZmqOtN3Tm1V3Q6gqttFpJbf/ebnuF/9nEGJyEhgJEDDhg0L8BglL3HvER6aspR2MVV4dEDeA91KgojQuVF1OjeqTtK+I/zv581M+HULny3bTseGVbmleywD2tQ5XiqYuiiJR6cup6a1LxhTbhUkOeTVkpmzLqofsAToDZwFfCMiPwLBuJJALLAf+EhErlfVD8/wfqjqm8Cb4KqV8nuIknY0M5t7JixGFUYP71RiYxDyE1MtgscuacWoPs2YsjCJd3/ayL0TFlOvSjg3nteYnSlpvPuTjV8wprwrSHJIAhr4vY7h96qhY0YAz6prwEgQkY1AS6ARsFFVkwFEZCpwHvAhsFNE6vpKDXWBXadxv4D33FdrWJq4nzeu6xQwU134qxQWwk3nNeaGbo2YvWYXY3/ayLNfrgGw9gVjTIGSwwKgmYjEAluBYcC1OY7ZAvQBfhSR2kALYAOuFNBNRCJw1Up9gGMtx9OBm4BnfX9+6rd9vIi8hGuQbgb8Wqin88ix+ZBuOrcRA9rW9TqcUwoKEvq2rk3f1rVZvT2FvYeP0r1pTa/DMsZ4LN/koKqZInIP8DWummisqq4UkTt8+8cATwPvichyXEJ4WFV3A7tFZAqwCNdAvRhfVRAuKUwWkVtxyWWI73orfb2aVvnOuftUPZUCTdK+Izz40VLa1K/MYwO9a2cojFZ1K3sdgjEmQNgI6SKUkZXNNf/9mXU7D/HZvT1oXLOS1yEZY8xJ2QjpEvL812tZvGU/o6/taInBGFOqWYtjEZm9Zidv/rCB67s1ZFC7el6HY4wxZ8SSQxHYtj+VByYvpXXdyjw+sGgm1DPGGC9ZcjhDGVnZ3DthMRmZ2bx2XSfCQ4O9DskYY86YtTmcoRdn/sbCzft4dXhHYq2dwRhTRljJ4Qx8t3YXY+asZ3jXhgxub+0Mxpiyw5JDIW0/kMoDk5bQsk4UT1xq7QzGmLLFkkMhZGZl88cJi0m3dgZjTBllbQ6F8O9Zv7Fg0z5eHtqBs6IjvQ7HGGOKnJUcTtOc35J5/fv1DI1rwOUdc80kbowxZYIlh9OwMyWNByYtoXmtKJ4cfLbX4RhjTLGx5FBAx9oZjhzN4rXrOlKxgrUzGGPKLmtzKKBXv13HLxv38uKQ9jStFeV1OMYYU6ys5FAAc9ft5j/fJXB15xiu6mxrKRtjyj5LDvnYlZLGfZMW0zQ6kr9dZu0MxpjywaqVTiErWxk1cQmH0jMZf1s3IirY22WMKR/s0+4UXv12HT9v2MO/rm5H89rWzmCMKT+sWukk5iXs5tXZ67iyY32GWDuDMaacseSQh+SD6YyatIQmNSvx9OVtEBGvQzLGmBJl1Uo5ZGUr909aQkpqBh/c2pVKYfYWGWPKnwKVHESkv4isFZEEEXkkj/1VRGSGiCwVkZUiMsK3vYWILPH7SRGR+3z7Jvlt3yQiS3zbG4tIqt++MUX5wPl57bsE5ibs5qnBZ9OyTuWSvLUxxgSMfL8Wi0gw8BpwEZAELBCR6aq6yu+wu4FVqnqpiEQDa0VknKquBTr4XWcrMA1AVYf63eNF4IDf9daraocze7TTN3/DHl6e9RuXd6jH0C4NSvr2xuR2ZC98chd0vA5aXep1NKYcKUjJoSuQoKobVPUoMBG4LMcxCkSJq5yPBPYCmTmO6YP70N/sv9F3zjXAhELEX2R2H0rnjxMW07hGJf5+RVtrZzCB4du/wW9fwqQb4Ne3vI7GlCMFSQ71gUS/10m+bf5GA62AbcByYJSqZuc4Zhh5J4CewE5VXee3LVZEFovIHBHpmVdQIjJSROJFJD45ObkAj3Fy2b52hv2pGYy+thOR1s5gAsHWhbDwPYi7BVoMgC8ehFlPgarXkZlyoCDJIa+v0Dn/dfYDlgD1cNVIo0XkeIW9iFQABgMf5XGt4ZyYNLYDDVW1I/AAMN7/WscDUH1TVeNUNS46OroAj3Fyb8xZz4/rdvPEpa1pXc/aGUwAyM6Cz/8EkbWg71NwzQfQ+WaY+xJ8cidkZXgdoSnjCpIckgD/CvgYXAnB3whgqjoJwEagpd/+AcAiVd3pf5KIhABXApOObVPVdFXd4/t9IbAeaF6wxzl9v27cy4sz13Jp+3pc27Vhcd3GmNOz6H3YthgufgbCK0NwCAx6GS78MyydAOOvgfSDXkdpyrCCJIcFQDMRifWVAIYB03McswXXpoCI1AZaABv89ucsHRzTF1ijqknHNohItK/xGhFpAjTLca0is+dQOvdOWETD6hH84wobz2ACxOHdrvqocU9oe/Xv20Xg/Idg8GjYMAfeGwgHd578OsacgXyTg6pmAvcAXwOrgcmqulJE7hCRO3yHPQ2cJyLLgW+Bh1V1N4CIROB6Ok3N4/J5tUP0ApaJyFJgCnCHqu49/UfL3/YDaYSFBDP62k5EhYcWxy2MOX2znoSjh+CSF1xCyKnTDTB8IuxeB+9cBLsTSjxEU/aJloHGrbi4OI2Pjy/UuZlZ2YQE20BxEyASf3Uf+Of9ES5++tTHJi2E8UPc79dOhpi44o/PlCkislBV8/yHU+4/FS0xmICRnQWfPwBR9eD8h/M/PqYz3PoNhEXBe4Ng7VfFH6NxMtIgdZ/XURQr+2Q0JlAseAd2LIf+/4CwyIKdU+MslyCiW8DE4bDw/eKN0UBaCrzdB0Z3gf2J+R9fSllyMCYQHNoFs/8OTS6A1pef3rmRteDmz6HJhTDjj/D9czYWorhkZcDkGyF5jSs9TBgORw97HVWxsORgTCD45q+QceTkjdD5CYuEaydB+2vh+3/AjFGQlXOSAnNGVGHGfbDhO7j0FRjyLuxaCdPugOycY35LP0sOxnht8zw3duG8e6Fms8JfJzgULn8dev7JjZOYdD0cPVJ0cZZ3PzwPSz6EXg9Bx+uh2UVw0dOwejrMedbr6IqcJQdjvJSVCZ8/CFUaQK8Hz/x6ItDnr64E8ttX8P6lcHjPmV+3vFs6Eb57BtoNgwsf+337uXdDh+thznOwIq/e+qWXJQdjvPTrm65qov8/oUKlortu19tg6Aewc4XrGrtvU9Fdu7zZMAc+vccNShz8nxOr/URg0EvQ4Bw3e+62Jd7FWcQsORjjlZTt8N0/oOlF0H32ygIAABxjSURBVHJQ0V+/1aVw46dwZA+8fVGZ+uAqMbtWuxlxazSFoR9CSIXcx4SEuX0RNWDitWVm1LolB2O88s1fIOsoDHiucI3QBdGwG9w6032AvTcQEr4tnvuURQd3wLghEBoO102GilVPfmxkLRg+wY19mHSd68lUyllyMMYLG3+A5R9Bj/vcWIXiFN3CjYWo1thN2Ld0YvHeryxIP+TeqyN73ejzqgWYlLNuO7jiv5C0wPUWK+XdiS05GFPSsjJcI3TVRtDj/pK5Z+W6MOILaHguTLsd5v671H94FZusTJgywg1IHPIe1DuNRSlbD3Yz5y6bCPNeLbYQS4IlB2NK2vzXYfdaGPAvCK1YcvcNrwLXfwxtrnaT+335kJuyw/xOFb78P1g30/X4an7x6V+j1//B2VfAN0+U6ilNLDkYU5IObHUjmJsPgBb9S/7+IWFw5Vtw7j2up9RHN5eJ+vEi89MrED8Wut8HXW4t3DVE4LLXXTXTx39wjdqlkCUHY0rS14+BZsEADwdNBQVBv2eg3z/cAK4Prijzk8gVyIqPYdYT0OYq6PPEmV2rQgQMm+D+nDDMtV2UMpYcjCkp62fDqk/cCOZqjb2Oxg3gunosbI2Hsf3hQFL+55RVm+e5aTAanuu+9QcVwUdjlfowdJzrsjz5xlK3tKslB2NKQmY6fPF/UL2JW6shULS5yrVDpGxzYyF2rvQ6opK3e52bQK9qQxg23nVdLSoNusDgV2HTj/BlAaZhDyCWHIwpCT+Phj0JMOD5ov3wKQqxveCWrwCFsQNg449eR1RyDiXDh1dBUAhcNwUiqhf9PdoPg+6jIP4d+PWtor9+MbHkYExx278F5jzvRiw36+t1NHmrfbYbCxFVBz68sszNE5Sno0dgwlA3Xfq1k6B6bPHdq88T0KyfKz1smFN89ylClhyMKW5fPep6sPT7p9eRnFrVBq4EUa8TTLkF5r/hdUTFJzsLpt4GWxfB1e8U/xKrQcFw1dtu1t3JN8Ke9cV7vyJgycGY4rTuG1jzmev7XrWB19HkL6I63PgJtBoEXz0CMx8vk2sV8PWf3d9L/2eh5cCSuWd4ZRg+0X1RmDAc0g6UzH0LqUDJQUT6i8haEUkQkUfy2F9FRGaIyFIRWSkiI3zbW4jIEr+fFBG5z7fvSRHZ6rfvEr/rPeq711oR6VdUD2tMicpIc43QNZq5cQWlRWhFGPI+dLkN5v0Hpo2EzKNeR1V0fn4dfnkDut0F3e4o2XtXj4Vr/gd717sxEAE8CDHf5CAiwcBrwACgNTBcRFrnOOxuYJWqtgcuAF4UkQqqulZVO6hqB6AzcASY5nfev4/tV9UvfPdrDQwDzgb6A6/7YjCmdPnpFdi3EQa+kPdsnoEsKBgued7VlS//CMZd7dZOLu1WTXdjTVoOgov/7k0Msb3c6Ph1M924igBVkJJDVyBBVTeo6lFgInBZjmMUiBIRASKBvUDONQr7AOtVdXM+97sMmKiq6aq6EUjwxWBM6bF3I8x9Cc6+0q0LXRqJQM8H4PIxsPkneHcAJMV7HVXhJS5w7Qz1O7tR4kEefufscuvvJbMl472L4xQKkhzqA4l+r5N82/yNBloB24DlwChVzVlROQyYkGPbPSKyTETGiki107ifMYHtq0dc98h+z3gdyZnrMNz15jm4A97uA+OugW2LvY7q9Ozd4HomRdV1z1IhwuuI3AJPsb3cDK5bfvE6mlwKkhzymmg+53SO/YAlQD2gAzBaRCofv4BIBWAw8JHfOW8AZ/mO3w68eBr3Q0RGiki8iMQnJycX4DGMKSFrv3RLdF7wCFSu53U0RaNpXxi11FUzJf0Kb14AE6+DHSu8jix/R/bCh1eDZruxDJVqeh2RExzq2nYq13drQOxPzP+cElSQ5JAE+HeziMGVEPyNAKaqkwBsBFr67R8ALFLV40skqepOVc3ylTDe4veqo4LcD1V9U1XjVDUuOjq6AI9hTAk4esTNdhrdCs4p4cbO4hYW6aqZRi1z01Jv/BHGdIfJN8GuNV5Hl7eMNNcz6ECS6ylUs6nXEZ0oororyWSmw8ThcPSw1xEdV5DksABoJiKxvhLAMGB6jmO24NoUEJHaQAtgg9/+4eSoUhKRun4vrwCOfQWZDgwTkTARiQWaAb8W7HGM8djcf7tBbwNfcN8My6LwynD+Q3DfUuj1kFtd7vVurvfN7gSvo/tddjZ8cgckzocr/+tWxQtE0S3gqndcKeyTOwOm63C+yUFVM4F7gK+B1cBkVV0pIneIyLGvRk8D54nIcuBb4GFV3Q0gIhHARUDOIZf/EpHlIrIMuBC433e/lcBkYBXwFXC3qgZufy9jjtmzHn56GdpeA417eB1N8atYDXr/Ge5b5la0W/M5vNYFpt3pGuS9NusJWDkNLvqbW18hkDW/GC5+GlZ9Cj/8y+toABAtA6tBxcXFaXx8Ke5FYUo/VTdHT9ICuCceomp7HVHJO5TskuOCtyE7Ezpc6xv8V4AlNovagrfh8z9Blz+4RXuKa43uoqQKn94NS8a5FehKIKGJyEJVzXN4uI2QNqYorJ4B67+FCx8rn4kBIDLa9c4atdR9KC+dCK92gs8ecIsclZS1X7nBh837Q//nSkdiABfnoH9DTFdX+tq+1NtwrORgzBk6ehhGd4WKVWHkHAgO8TqiwHBgK/z4Iiz6H0gQxI1wa2ZH1Sm+e25bDO9eAjWbuzWzK1QqvnsVl4M74a3egMJt3xXrlw0rORhTnH54HlKSYOCLlhj8VakPg16CexdC+6FuuupXOrh5jQ7vLvr77dsM44dCRE24dnLpTAzgksHw8W51vknXebaMqyUHY85E8m8wbzR0uC5we8N4rVojGPwfuDfe1aPPfx1ebgezniy65TNT98G4Ie6D9LqPSn/VXt32cPkbrg3rs/tce0QJs+RgTGGpwhcPutG2fZ/yOprAV70JXPEG3P0rtBgAc192SWL2M5C6v/DXzUyHSTe4UdDDxkGtlvmfUxqcfTlc8CgsneCm2ShhlhyMKayV02DjHOj9F9cYawqmZjO3hsJdP0PT3q7r5svtYM6/Tn9yP1WYfq9bhvPy1yG2Z/HE7JVeD0Hry+Gbv8JvM0v01pYcjCmM9INuds+67SHuFq+jKZ1qtXLTV98x140L+e4ZeKUd/PgSpB8q2DW+ewaWTYLej0O7a4o3Xi8EBbnqpTpt3QJMJTgS3ZKDMYUx5zk4uB0GvuTt7J5lQZ22rgF25PcQ0wW+fQpeae+qUo4eOfl5iz5wnQE63Qg9HyypaEtehQgYPsGtszFhWNG10+TDkoMJPKtnwM+vBcYo27zsWu2W0Ox0Y/EvL1me1OvoGpNvneUSxszH4dUOMH9M7h47CbPcbKZn9XEJurSMZSisKjGuPSVlq1tmNCuj2G9p4xxMYFk5DT4awfGJeGu3dUtWtroUarX2/kNAFd4bBLtWwj0LoVINb+MpyzbPg+/+4doToupBrz9Bxxth91oY2x+qNYYRX7q5nsqLJRPcfFFd/uC6Tp+hU41zsORgAkfCt66fekwcDHrZfTtc8xlsmQ+o6+3SchC0GuwWbAnyoOC7bLJbMGbQy25Qlyl+G+a4toXEX6BKA/etOSgY/jCr7EyJfjpm/gXmveqSQ5c/nNGlLDmYwJe4AP432CWAmz93o42PObgT1n4Oqz9zvYOyM92iLS0HuhJFo+4lMwNq2gEY3cUV8W+d5U1yKq9UYf1slyT2JMDNX0CdNl5H5Y3sLDcNecIsuPETt2BQIVlyMIFt12q3BGV4Vbjl61MPYErd79beXT3dlTQyjrjzWlziqp/O6u0a7orDl4/AL2PgttlQv1Px3MOcmqorOZS2NbmLWloKvN0XDu9y/x6rNynUZSw5mMC1bzOM7ef+09/yFVSPLfi5R4+4b5OrZ8BvX7pv9qERbtWyVoPdNMjhVYomzh3L4b+9oPMINyWEMV7bs94t21qvE9yQc0WEgjlVcrCJYIx3Du2CDy533/5HfHl6iQFcF79Wg9xPVoZruFz9mWunWD0dgkKhyfmu6qnFwMIPVMvOhs8f9K1f8HjhrmFMUatxllv2tJimRLeSg/FG2gHX62f3OrjxU2h4TtFdOzvbzUmzZoYrVezbBAg0PNclilaDTu8/1JLxboWuy16DjtcXXZzGeMyqlUxgyUh1C+Mk/gLDJ0GzvsV3L1XYudIliTWfwU7farR120PLS12yiG5x8i6yqfvgP3GuTveWr60R2pQpVq1kAkdWppsGYPM8uOrt4k0M4D7067RxPxc+6upp13zmqp+++7v7qdHUV6K41NXf+ieK2c9A6l4YOM0SgylXrORgSk52tlsGcel4t3Rj19u8jSdlu6+L7AzY+CNoFlSu7xtLcalr03i7L3S5DS4JjHV9jSlKVq1kvKfqFnmZ/xpc8Bhc8LDXEZ3oyF747evfl/vM9E3XUKkW3LPgxHEXxpQRVq1kvPfjiy4xdL0dzn/I62hyi6gOHYa7n6OH3QCjdTPd4jSWGEw5VKBKVBHpLyJrRSRBRB7JY38VEZkhIktFZKWIjPBtbyEiS/x+UkTkPt++50VkjYgsE5FpIlLVt72xiKT6nTOmKB/YeCB+LMx+GtpeA/2f9X5+pPxUqAStL3O9k5oWc5uIMQEq3+QgIsHAa8AAoDUwXERa5zjsbmCVqrYHLgBeFJEKqrpWVTuoagegM3AEmOY75xugjaq2A34DHvW73vpj56nqHWfwfMZrK6fBZw9As35uMRZr1DWmVCjI/9SuQIKqblDVo8BE4LIcxygQJSICRAJ7gcwcx/TBfehvBlDVmap67Jj5QEwhn8EEqoRv4ePb3NrKQ94rmfmPjDFFoiDJoT6Q6Pc6ybfN32igFbANWA6MUtXsHMcMAyac5B63AF/6vY4VkcUiMkdE8lz3T0RGiki8iMQnJycX4DFMiUpcAJOud2MIhk90PX+MMaVGQZJDXhXEObs49QOWAPWADsBoETk+ybqIVAAGAx/lurjIn3GljHG+TduBhqraEXgAGO9/reMBqL6pqnGqGhcdbev3BpRdq2H8EIisDddPtQZdY0qhgiSHJKCB3+sYXAnB3whgqjoJwEagpd/+AcAiVd3pf5KI3AQMAq5TX59aVU1X1T2+3xcC64HmBX8k46l9m+GDKyA4DG6YduoZVo0xAasgyWEB0ExEYn0lgGHA9BzHbMG1KSAitYEWwAa//cPJUaUkIv2Bh4HBqnrEb3u0rxEcEWkCNMtxLROoDu1yiSHjiJsl8nQn0jPGBIx8xzmoaqaI3AN8DQQDY1V1pYjc4ds/BngaeE9EluOqoR5W1d0AIhIBXATcnuPSo4Ew4BvXjs18X8+kXsDfRCQTyALuUNWSWVHbFF7aATdfUso2N5Fe7bO9jsgYcwZshLQ5cyU5kZ4xpsjYCGlTfEp6Ij1jTImw5GAKLzsbpt8La79wE+m1vdrriIwxRcSGq5rCUYWZj7sZVi94zPsZVo0xRcqSgymcYxPpnXNHYE6kZ4w5I5YczOnzn0iv3z8DfyI9Y8xps+RgTo9NpGdMuWD/s03BrZ9tE+kZU05YcjAFkxQPE20iPWPKC0sOJn+7VsO4qyGylk2kZ0w5YcnBnJr/RHo3fmIT6RlTTtggOHNy/hPpjfgSqjX2OiJjTAmx5GDyZhPpGVOuWXIwuWWkwoThsGuVm0iv4TleR2SMKWGWHMyJbCI9YwyWHIw/m0jPGONjycE4+zbD14/Bms/gwj/bRHrGlHOWHMq79EMw998w7z8gQdD3Keg+yuuojDEes+RQXmVnw7KJMOspOLTDTaLX90moUt/ryIwxAcCSQ3m05Rf46hHYtgjqd4ahH0KDLl5HZYwJIJYcypMDSfDNE7BiCkTVhSvehLZDbGZVY0wuBfpUEJH+IrJWRBJE5JE89lcRkRkislREVorICN/2FiKyxO8nRUTu8+2rLiLfiMg635/V/K73qO9ea0WkX1E9bLl19DB89w/4T5xrcO71ENy7ENoPtcRgjMlTviUHEQkGXgMuApKABSIyXVVX+R12N7BKVS8VkWhgrYiMU9W1QAe/62wFpvnOeQT4VlWf9SWcR4CHRaQ1MAw4G6gHzBKR5qqaVRQPXK6owvKPXGnh4DY4+0q46Cmo2tDryIwxAa4g1UpdgQRV3QAgIhOBywD/5KBAlIgIEAnsBTJzXKcPsF5VN/teXwZc4Pv9feB74GHf9omqmg5sFJEEXww/n9aTlXdJC+GrhyFpAdTtAFePhUbneh2VMaaUKEhyqA8k+r1OAnLOpzAamA5sA6KAoaqaneOYYcAEv9e1VXU7gKpuF5Fafvebn+N+ubrQiMhIYCRAw4b2Tfi4lG2uB9KyiRBZGy57HdoPt+ojY8xpKUhyyGuBYM3xuh+wBOgNnAV8IyI/qmoKgIhUAAYDjxbR/VDVN4E3AeLi4nLtL3cyUt1Yhbn/huws6PEA9HwAwqK8jswYUwoVJDkkAQ38XsfgSgj+RgDPqqoCCSKyEWgJ/OrbPwBYpKo7/c7ZKSJ1faWGusCu07ifOUYVVk517QoHEqHVYLj4aZte2xhzRgpS17AAaCYisb4SwDBcFZK/Lbg2BUSkNtAC2OC3fzgnVinhu8ZNvt9vAj712z5MRMJEJBZoxu9JxvjbugjG9ncT5YVXhZs/h6EfWGIwxpyxfEsOqpopIvcAXwPBwFhVXSkid/j2jwGeBt4TkeW4aqGHVXU3gIhE4Ho63Z7j0s8Ck0XkVlxyGeK73koRmYxr8M4E7raeSjkc3AHf/g2WjIdKNeHSV6Hj9RAU7HVkxpgyQlxNUOkWFxen8fHxXodR/DLS4OfR8ONLkHUUut0Jvf4Pwit7HZkxphQSkYWqGpfXPhshXRqowqpP4Zu/wP4t0GKga1eocZbXkRljyihLDoFu+1L46jHYPBdqne2W7GxygddRGWPKOEsOgerQLpj9NCz6ACKqw8CXoNNNEGx/ZcaY4mefNIEmMx3mvwE/vACZqdDtLjj/IahY1evIjDHliCWHQKEKaz6HmY/Dvo3QvD9c/AzUbOp1ZMaYcsiSQyDYnQCf3w8bf4DolnD9VGjax+uojDHlmCUHL6lC/FhXWgiuAAOeh7hbrF3BGOM5+xTyyqFd8Ok9sO5raHIhXP4GVK7rdVTGGANYcvDGmi9g+r2QfhD6PwddR9qsqcaYgGLJoSSlH4KvH4NF70OdtnDlZ1CrlddRGWNMLpYcSkpSPEy9DfZuhO73wYWPQUiY11EZY0yeLDkUt6xM+PEFmPMvqFwPbv4MGvfwOipjjDklSw7Fac96mDoStsZDu2Fwyb8gvIrXURljTL4sORQHVdeu8NWjEBzq1m9uc5XXURljTIFZcihqh5Jhxh9h7RcQe77rolol1xLYxhgT0Cw5FKXfvoZP74a0FOj3DzjnTuuiaowplSw5FIWjh90o5/ixULsN3Dgdarf2OipjjCk0Sw5nautC1+i8Zz2cdy/0/ot1UTXGlHqWHAorKxPmvgRznoPIOnDTdIjt5XVUxhhTJCw5FMbeDTD1dkj6FdoOgUtesPUWjDFlSoFaS0Wkv4isFZEEEXkkj/1VRGSGiCwVkZUiMsJvX1URmSIia0RktYic69s+SUSW+H42icgS3/bGIpLqt29MUT3sGVN1K7ON6QnJa+Gqd+Cqty0xGGPKnHxLDiISDLwGXAQkAQtEZLqqrvI77G5glapeKiLRwFoRGaeqR4FXgK9U9WoRqQBEAKjqUL97vAgc8LveelXtcKYPV6QO73FdVNd8Bo17whVjoEqM11EZY0yxKEi1UlcgQVU3AIjIROAywD85KBAlIgJEAnuBTBGpDPQCbgbwJYuj/hf3nXMN0PuMnqQ4rfvGdVFN3QcX/x263W1dVI0xZVpBPuHqA4l+r5N82/yNBloB24DlwChVzQaaAMnAuyKyWETeFpFKOc7tCexU1XV+22J9x88RkZ55BSUiI0UkXkTik5OTC/AYhXD0CHz+IIy7GiJqwG3fuR5JlhiMMWVcQT7lJI9tmuN1P2AJUA/oAIz2lRpCgE7AG6raETgM5GyzGA5M8Hu9HWjoO/4BYLzvWicGoPqmqsapalx0dHQBHuM0bVsMb54PC96Cc+9xiaFOm6K/jzHGBKCCJIckoIHf6xhcCcHfCGCqOgnARqCl79wkVf3Fd9wUXLIAQERCgCuBSce2qWq6qu7x/b4QWA80P52HOiPZWfDDC/B2X7f+wo2fQr9nIDS8xEIwxhivFSQ5LACaiUisr0F5GDA9xzFbgD4AIlIbaAFsUNUdQKKItPAd14cT2yr6AmtUNenYBhGJ9jWCIyJNgGbAhtN+ssLYtwnevQRmPw2tBsNd86DJBSVya2OMCST5NkiraqaI3AN8DQQDY1V1pYjc4ds/BngaeE9EluOqoR5W1d2+S9wLjPMllg24UsYxwzixSglcA/bfRCQTyALuUNW9hX7CglCFpRPgi4dABK58y41fkLxq1IwxpuwT1ZzNB6VPXFycxsfHF+7kI3thxihYPR0a9YAr3oCqDYs2QGOMCUAislBV4/LaV75HSG9bDOOHwZE9cNHfXMNzULDXURljjOfKd3Ko2ghqtXKJoW47r6MxxpiAUb6TQ0R1uPETr6MwxpiAY6O5jDHG5GLJwRhjTC6WHIwxxuRiycEYY0wulhyMMcbkYsnBGGNMLpYcjDHG5GLJwRhjTC5lYm4lEUkGNp/BJWoCu/M9qnyw9+JE9n78zt6LE5WF96ORqua5IE6ZSA5nSkTiTzb5VHlj78WJ7P34nb0XJyrr74dVKxljjMnFkoMxxphcLDk4b3odQACx9+JE9n78zt6LE5Xp98PaHIwxxuRiJQdjjDG5WHIwxhiTS7lODiLSX0TWikiCiDzidTxeEpEGIvKdiKwWkZUiMsrrmLwmIsEislhEPvM6Fq+JSFURmSIia3z/Rs71OiYvicj9vv8nK0RkgoiEex1TUSu3yUFEgoHXgAFAa2C4iLT2NipPZQJ/UtVWQDfg7nL+fgCMAlZ7HUSAeAX4SlVbAu0px++LiNQH/gjEqWobIBgY5m1URa/cJgegK5CgqhtU9SgwEbjM45g8o6rbVXWR7/eDuP/89b2NyjsiEgMMBN72OhaviUhloBfwDoCqHlXV/d5G5bkQoKKIhAARwDaP4yly5Tk51AcS/V4nUY4/DP2JSGOgI/CLt5F46mXgISDb60ACQBMgGXjXV832tohU8joor6jqVuAFYAuwHTigqjO9jarolefkIHlsK/f9ekUkEvgYuE9VU7yOxwsiMgjYpaoLvY4lQIQAnYA3VLUjcBgot210IlINV8sQC9QDKonI9d5GVfTKc3JIAhr4vY6hDBYNT4eIhOISwzhVnep1PB7qDgwWkU246sbeIvKhtyF5KglIUtVjJckpuGRRXvUFNqpqsqpmAFOB8zyOqciV5+SwAGgmIrEiUgHXoDTd45g8IyKCq1NeraoveR2Pl1T1UVWNUdXGuH8Xs1W1zH0zLChV3QEkikgL36Y+wCoPQ/LaFqCbiET4/t/0oQw20Id4HYBXVDVTRO4Bvsb1Nhirqis9DstL3YEbgOUissS37TFV/cLDmEzguBcY5/sitQEY4XE8nlHVX0RkCrAI18tvMWVwKg2bPsMYY0wu5blayRhjzElYcjDGGJOLJQdjjDG5WHIwxhiTiyUHY4wxuVhyMMYYk4slB2OMMbn8P+5We/BcPyVFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Loss curve\n",
    "plt.plot(train_loss)\n",
    "plt.plot(dev_loss)\n",
    "plt.title('Loss')\n",
    "plt.legend(['train', 'dev'])\n",
    "plt.savefig('loss.png')\n",
    "plt.show()\n",
    "\n",
    "# Accuracy curve\n",
    "plt.plot(train_acc)\n",
    "plt.plot(dev_acc)\n",
    "plt.title('Accuracy')\n",
    "plt.legend(['train', 'dev'])\n",
    "plt.savefig('acc.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HzIcYAfvkUZ_"
   },
   "source": [
    "### Predicting testing labels\n",
    "\n",
    "Predictions are saved to *output_logistic.csv*.\n",
    "\n",
    "預測測試集的資料標籤並且存在 *output_logistic.csv* 中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 187
    },
    "colab_type": "code",
    "id": "ZEAKhugPkUyH",
    "outputId": "d18595c8-7851-4f76-e04e-51accf15bc53"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Not in universe -4.031960278019251\n",
      " Spouse of householder -1.625403958705141\n",
      " Other Rel <18 never married RP of subfamily -1.4195759775765404\n",
      " Child 18+ ever marr Not in a subfamily -1.2958572076664745\n",
      " Unemployed full-time 1.1712558285885912\n",
      " Other Rel <18 ever marr RP of subfamily -1.167791807296237\n",
      " Italy -1.093458143800618\n",
      " Vietnam -1.0630365633146415\n",
      "num persons worked for employer 0.938992277356651\n",
      " 1 0.8226614922117185\n"
     ]
    }
   ],
   "source": [
    "# Predict testing labels\n",
    "predictions = _predict(X_test, w, b)\n",
    "with open(output_fpath.format('logistic'), 'w') as f:\n",
    "    f.write('id,label\\n')\n",
    "    for i, label in  enumerate(predictions):\n",
    "        f.write('{},{}\\n'.format(i, label))\n",
    "\n",
    "# Print out the most significant weights\n",
    "ind = np.argsort(np.abs(w))[::-1]\n",
    "with open(X_test_fpath) as f:\n",
    "    content = f.readline().strip('\\n').split(',')\n",
    "features = np.array(content)\n",
    "for i in ind[0:10]:\n",
    "    print(features[i], w[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1C6sqhUbLMGe"
   },
   "source": [
    "# Porbabilistic generative model\n",
    "\n",
    "In this section we will discuss a generative approach to binary classification. Again, we will not go through the formulation detailedly. Please find [Prof. Lee's lecture](https://www.youtube.com/watch?v=fZAZUYEeIMg) if you are interested in it.\n",
    "\n",
    "接者我們將實作基於 generative model 的二元分類器，理論細節請參考[李宏毅老師的教學影片](https://www.youtube.com/watch?v=fZAZUYEeIMg)。\n",
    "\n",
    "### Preparing Data\n",
    "\n",
    "Training and testing data is loaded and normalized as in logistic regression. However, since LDA is a deterministic algorithm, there is no need to build a development set.\n",
    "\n",
    "訓練集與測試集的處理方法跟 logistic regression 一模一樣，然而因為 generative model 有可解析的最佳解，因此不必使用到 development set。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "czWXO7qML8DU"
   },
   "outputs": [],
   "source": [
    "# Parse csv files to numpy array\n",
    "with open(X_train_fpath) as f:\n",
    "    next(f)\n",
    "    X_train = np.array([line.strip('\\n').split(',')[1:] for line in f], dtype = float)\n",
    "with open(Y_train_fpath) as f:\n",
    "    next(f)\n",
    "    Y_train = np.array([line.strip('\\n').split(',')[1] for line in f], dtype = float)\n",
    "with open(X_test_fpath) as f:\n",
    "    next(f)\n",
    "    X_test = np.array([line.strip('\\n').split(',')[1:] for line in f], dtype = float)\n",
    "\n",
    "# Normalize training and testing data\n",
    "X_train, X_mean, X_std = _normalize(X_train, train = True)\n",
    "X_test, _, _= _normalize(X_test, train = False, specified_column = None, X_mean = X_mean, X_std = X_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L8NGKl-EPvok"
   },
   "source": [
    "### Mean and Covariance\n",
    "\n",
    "In generative model, in-class mean and covariance are needed.\n",
    "\n",
    "在 generative model 中，我們需要分別計算兩個類別內的資料平均與共變異。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "iQrzXXKUPwHT"
   },
   "outputs": [],
   "source": [
    "# Compute in-class mean\n",
    "X_train_0 = np.array([x for x, y in zip(X_train, Y_train) if y == 0])\n",
    "X_train_1 = np.array([x for x, y in zip(X_train, Y_train) if y == 1])\n",
    "\n",
    "mean_0 = np.mean(X_train_0, axis = 0)\n",
    "mean_1 = np.mean(X_train_1, axis = 0)  \n",
    "\n",
    "# Compute in-class covariance\n",
    "cov_0 = np.zeros((data_dim, data_dim))\n",
    "cov_1 = np.zeros((data_dim, data_dim))\n",
    "\n",
    "for x in X_train_0:\n",
    "    cov_0 += np.dot(np.transpose([x - mean_0]), [x - mean_0]) / X_train_0.shape[0]\n",
    "for x in X_train_1:\n",
    "    cov_1 += np.dot(np.transpose([x - mean_1]), [x - mean_1]) / X_train_1.shape[0]\n",
    "\n",
    "# Shared covariance is taken as a weighted average of individual in-class covariance.\n",
    "cov = (cov_0 * X_train_0.shape[0] + cov_1 * X_train_1.shape[0]) / (X_train_0.shape[0] + X_train_1.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kifW1pFxXXA5"
   },
   "source": [
    "### Computing weights and bias\n",
    "\n",
    "Directly compute weights and bias from in-class mean and shared variance. [Prof. Lee's lecture slides](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Classification%20(v3).pdf)(p.33) gives a concise explanation.\n",
    "\n",
    "權重矩陣與偏差向量可以直接被計算出來，算法可以參考[李宏毅老師教學投影片](http://speech.ee.ntu.edu.tw/~tlkagk/courses/ML_2016/Lecture/Classification%20(v3).pdf)第 33 頁。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "UghOxYrUXXPU",
    "outputId": "72f8d4b8-0019-414e-f7d0-9dfc4796188d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training accuracy: 0.8637017104099086\n"
     ]
    }
   ],
   "source": [
    "# Compute inverse of covariance matrix.\n",
    "# Since covariance matrix may be nearly singular, np.linalg.inv() may give a large numerical error.\n",
    "# Via SVD decomposition, one can get matrix inverse efficiently and accurately.\n",
    "u, s, v = np.linalg.svd(cov, full_matrices=False)\n",
    "inv = np.matmul(v.T * 1 / s, u.T)\n",
    "\n",
    "# Directly compute weights and bias\n",
    "w = np.dot(inv, mean_0 - mean_1)\n",
    "b =  (-0.5) * np.dot(mean_0, np.dot(inv, mean_0)) + 0.5 * np.dot(mean_1, np.dot(inv, mean_1))\\\n",
    "    + np.log(float(X_train_0.shape[0]) / X_train_1.shape[0]) \n",
    "\n",
    "# Compute accuracy on training set\n",
    "Y_train_pred = 1 - _predict(X_train, w, b)\n",
    "print('Training accuracy: {}'.format(_accuracy(Y_train_pred, Y_train)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RDKWzBy0bi3c"
   },
   "source": [
    "### Predicting testing labels\n",
    "\n",
    "Predictions are saved to *output_generative.csv*.\n",
    "\n",
    "預測測試集的資料標籤並且存在 *output_generative.csv* 中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 187
    },
    "colab_type": "code",
    "id": "T3QjToT_Sq9J",
    "outputId": "0e22d03e-f1d4-49c0-f5db-6aaa4ad612f4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Retail trade 8.6865234375\n",
      " 34 -6.5498046875\n",
      " 37 -5.9912109375\n",
      " Other service -5.671875\n",
      " 19 5.361328125\n",
      " Transportation -5.3017578125\n",
      " 3 -5.29296875\n",
      " 46 4.40625\n",
      " Not in universe or children 4.400390625\n",
      " 9 -4.357421875\n"
     ]
    }
   ],
   "source": [
    "# Predict testing labels\n",
    "predictions = 1 - _predict(X_test, w, b)\n",
    "with open(output_fpath.format('generative'), 'w') as f:\n",
    "    f.write('id,label\\n')\n",
    "    for i, label in  enumerate(predictions):\n",
    "        f.write('{},{}\\n'.format(i, label))\n",
    "\n",
    "# Print out the most significant weights\n",
    "ind = np.argsort(np.abs(w))[::-1]\n",
    "with open(X_test_fpath) as f:\n",
    "    content = f.readline().strip('\\n').split(',')\n",
    "features = np.array(content)\n",
    "for i in ind[0:10]:\n",
    "    print(features[i], w[i])\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "hw2_classification.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3.6.8 64-bit ('ml': conda)",
   "language": "python",
   "name": "python36864bitmlconda4727b916de1c4a8f8265036fed6d2bb8"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
